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SECTION  I 
GRID  GENERATION 


1.  Background 

In  generating  a  grid  about  aircraft-like  configurations  several  constraints  should  be 
kept  in  mind.  To  begin,  the  grids  must  be  smoothly  varying  so  as  to  maintain  solu¬ 
tion  accuracy.  The  grids  should  also  be  body  conforming  to  enhance  solution  accuracy, 
simplify  the  implementation  of  boundary  conditions,  and  to  minimize  programming  com¬ 
plexity.  Finally,  in  order  to  use  approximately  factored  implicit  schemes  and  to  maintain 
computational  compatibility  with  vectorized  machines,  the  grids  should  be  well-ordered. 

For  a  simple  wing-body  combination  one  can  envision  a  single  mapping  procedure  as 
illustrated  in  Figure  (1).  This  warped  spherical  coordinate  system  is  well-ordered,  maps 
the  body  onto  a  single  coordinate  surface,  and  if  properly  generated,  it  can  be  sufficiently 
smooth.  Moreover,  the  axis  singularity  is  not  a  problem  for  Euler  and  thin  layer  Navier- 
Stokes  equations  that  are  transformed  in  general  coordinates  and  strong  conservation  law 
form. 

A  single  mapping  such  as  that  illustrated  in  Figure  (1)  is,  of  course,  inadequate  for 
complex  aerodynamic  configurations  which  can  include  engine  nacelles,  stores,  etc.  For 
these  cases,  subgrids,  which  are  either  embedded-to  or  overset-on  the  main  wing-body 
conforming  grid,  are  envisioned.  Some  possible  grid  configurations  are  illustrated  in  Figures 
(2)  and  (3)  for  two  dimensional  cross  sections. 

Embedding  or  oversetting  meshes  to  account  for  complex  configurations  requires  an 
immense  effort  in  interfacing  grids  and  numerical  algorithm  developments.  Consequently, 
this  university  contract  was  restricted  to  the  more  fundamental  task  of  generating  only  the 
main  body-conforming  grid  using  elliptic  partial  differential  grid  generation  equations  with 
clustering  terms.  We  have,  however,  begun  development  of  a  three  dimensional  hyperbolic 
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Figure  1.  Well-Ordered  Warped  Spherical  Grid  Mapping 


Figure  2.  Example  of  Mesh  Embedding  in  Two  Dimensions 


Figure  3.  Example  of  Mesh  Overlapping  in  Two  Dimensions 


grid  generation  procedure,  and,  in  two  dimensions,  we  have  indeed  studied  overset  grids. 

2.  Grid  Generation 

As  noted  earlier,  we  wish  to  build  a  single  main  grid  that  is  smoothly  continuous, 
maps  the  body  onto  a  constant  coordinate  surface,  and  is  well-ordered.  The  resulting 
mesh  must  also  use  surface  as  well  as  field  grid  points  efficiently.  Clustering  of  grid  points 
near  a  body  to  resolve  viscous  terms  should  not  cause  unwanted  extraneous  points  in 
the  free  stream.  These  considerations  lead  to  the  use  of  warped  spherical  and  cylindrical 
coordinate  systems. 

To  construct  this  main  spherical  grid  we  have  developed  grid  generation  codes  using 
both  elliptic  and  hyperbolic  partial  differential  equations.  The  elliptic  grid  generation 
procedure  has  been  written  up  in  the  form  of  a  technical  paper  and  was  presented  by  Mr. 
Reese  Sorenson  at  an  ASME  specialist  meeting.  This  paper,  which  describes  the  method 
and  shows  results,  is  reproduced  in  Appendix  A. 

A  paper  describing  our  work  in  three  dimensional  hyperbolic  grid  generation  procedure 
is  in  preparation.  Details  of  this  method  are  described  in  an  abbreviated  write-up  as 
Appendix  B. 

In  general  the  elliptic  grid  generator  is  more  powerful  as  it  can  generate  a  smooth  grid 
interior  to  user  specified  inner  and  outer  boundaries.  The  hyperbolic  method  generates 
a  grid  exterior  to  a  user  specified  inner  boundary.  The  outer  boundary  location  is  not 
known  in  advance  (although  one  could  iterate  on  its  location,  but  not  its  placement  of 
points).  However,  the  hyperbolic  solver  can  generate  an  orthogonal  grid  and  it  is  a  much 
more  efficient  code  to  use  than  the  elliptic  solver.  For  cases  in  which  it  is  suited,  it  also 
requires  less  user  input.  The  elliptic  solver,  however,  is  more  reliable. 

S.  Overset  Grids 

A  paper  describing  a  chimera  or  overset  grid  scheme  in  two  dimensions  is  attached  as 


Appendix  C.  This  paper  was  presented  at  the  same  ASME  meeting  as  the  three  dimensional 
elliptic  solver.  The  various  results,  which  are  obtained  with  linear  incompressible  flow 
equations,  verify  the  feasibility  of  this  approach,  but  considerably  more  research  is  needed. 


SECTION  n 
ZONAL  METHODS 


1.  Background 

Just  as  one  wants  to  limit  the  number  of  grid  points  to  only  those  action  regions 
in  which  they  are  needed,  one  would  also  like  to  limit  the  complexity  of  the  governing 
partial  differential  equations  to  appropriate  action  zones.  The  Navier-Stokes  equations, 
for  example,  are  not  needed  to  describe  the  flow  in  the  inviscid  far  field.  Indeed,  even 
nonlinear  full  potential  theory  is  not  needed  there,  as  a  linear  potential  approximation  is 
quite  adequate. 

In  principle,  considerable  computational  work  can  be  saved  by  using  just  that  s' 
plified  governing  equation  set  that  suffices  for  a  given  region  of  the  flow  field.  Compv 
storage  can  also  be  reduced.  However,  such  a  zonal  method  is  not  without  pitfalls.  .  >r 
one,  programming  and  matching  together  several  numerical  methods  increases  tue  overall 
complexity  of  the  computer  code  and  its  data  base.  This  is  especially  undersirable  on 
parallel  processors.  Moreover,  unless  one  is  very  careful  in  matching  the  different  numer¬ 
ical  algorithms  and  governing  equations,  overall  stability  and  iterative  convergence  can 
be  impaired.  Numerical  stability  could  conceivably  decline  so  much  that  the  simplified 
zonal  scheme  could  be  more  inefficient  than  a  straight  Navier-Stokes  algorithm  with  a  well 
stretched  grid. 

Such  negative  arguments  may  be  offset  by  additional  advantages  for  zonal  schemes. 
For  one,  steady  state  numerical  algorithms  for  potential  flow  are  far  advanced  over  those  for 
Euler  or  Navier-Stokes  equations.  This  is  because  the  simpler  scalar  potential  equation 
is  much  easier  to  optimize  for  steady  state  convergence.  A  well  matched  zonal  method 
may  therefore  have  better  convergence  than,  say,  a  Navier-Stokes  scheme  alone,  simply 
because  the  outer  flow  region  can  be  converged  at  a  much  faster  rate.  Another,  albeit 
intangible  aspect  of  the  zonal  method  is  that  one  is  forced  to  keep  the  numerical  algorithms 


compatible,  and  if  the  code  is  properly  designed,  either  algorithm  could  be  readily  used  in 
a  stand  alone  mode. 

Arguments  for  and  against  zonal  methods  will  be  summarized  later.  Clearly,  a  zonal 
scheme  has  some  reduced  computer  storage  over  a  Navier-Stokes  alone  algorithm,  and  it 
may  have  other  significant  computational  benefits.  Consequently,  we  have  coded  and  tested 
a  zonal  algorithm.  In  our  tests  we  have  essentially  combined  a  transonic  full  potential  code 
with  an  Euler  or  thin  layer  Navier-Stokes  code.  The  zonal  code  was  tested  in  first  two  and 
then  three  dimensions. 

2.  The  Zonal  Codes 

The  three  dimensional  zonal  code  which  we  have  written  matches  an  unsteady  con¬ 
servative  full  potential  code  with  an  Euler  or  thin  layer  Navier-Stokes  code.  Details  of  the 
full  potential  code,  partially  developed  for  this  application,  were  presented  as  an  AIAA 
paper  attached  as  Appendix  D.  The  zonal  code  itself  has  not  been  published,  but  Mr. 
Jack  Striegberger  is  using  the  code  in  his  Ph.D.  thesis  project,  so  it  will  ultimately  be  fully 
documented.  Mr.  Timothy  Barth  has  transferred  the  code  to  AFFDL. 

A  two  dimensional  zonal  code  has  also  been  written  which  combines  a  flux  vector 
split  Euler  code  with  a  potential  and  vector  potential  code.  Using  this  dual  potential 
combination,  the  outer  flow  is  able  to  correctly  convect  entropy.  This  has  lead  to  a  very 
versatile  zonal  code  in  which  we  are  able  to  match  the  equation  zones  in  a  much  more 
elegant  (and  simpler)  way.  Details  of  this  method  have  also  been  presented  at  an  AIAA 
technical  meeting,  and  this  paper  is  presented  as  Appendix  E. 

S.  Reflections 

The  zonal  codes  that  we  have  written  work  and  save  both  computer  time  and  some 
storage.  They  are,  however,  more  complex  than  a  single  Euler  or  Navier-Stokes  code. 
Moreover,  the  zonal  codes  are  not  as  readily  generalized  to  a  new  problem  because  some¬ 
where  in  the  flow  field  they  use  simplifying  assumptions.  In  fact,  a  zonal  code  makes  a 
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trade  between  computer  time  and  engineering  time.  A  zonal  code  that  gives  the  same 
result  as  a  Navier-Stokes  code  will  be  cheaper  to  run  on  the  computer,  but  it  will  require 
more  engineering  development  time.  So  the  zonal  code  is  ideal  for  optimizing  a  given  con* 
figuration  because,  once  the  code  is  set  up,  it  will  run  more  efficiently.  If  one  is  continually 
changing  the  layout  of  the  configuration  and  the  type  of  flow  field  being  solved,  then  a 
nonzonal  general  code  may  be  more  economical  because  the  engineering  time  for  the  first 
solution  will  be  less. 

The  area  in  which  more  research  is  needed  with  zonal  codes  is  in  how  to  interface 
the  zones  more  tightly  with  a  minimum  of  bookkeeping.  The  two  dimensional  zonal  code 
described  in  Appendix  E  presents  one  approach  that  I  believe  we  can  and  should  generalize. 

In  the  zonal  method  described  in  Appendix  E  we  solve  (using  u  =  <Pt  +  rpr,v  =  <f>,-rpt) 


(pu)j  +  (Pv)»  =  0  (la) 

V,  U,  =  uj  (16) 

us,  +  vs,  =0  (lc) 

p  =  p{u2  +  v2,s)  (Id) 


throughout  the  entire  flow.  These  equations  can  convect  but  not  produce  entropy.  In  an 
entropy  producing  zone  (e.g.,  around  a  shock,  see  Figure  (lb)  of  Appendix  E)  the  inviscid 
conservation  law  equations  of  mass,  momentum,  and  energy  are  solved.  From  this  solution 
we  obtain  the  value  of  w  that  feeds  into  Equation  (lb)  by  forming  w  =  (u,  -  v,)  directly 
from  the  conservation  law  solutions.  Elsewhere  u  is  evaluated  from 

-u)  =  (7Af2)_1(us,  -  us,) 
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so  that  Equation  (lb)  is  the  Crocco  equation.  Likewise  in  the  shock  zone  a  is  overloaded 
directly  from  the  conservation  laws  in  place  of  Equation  (lc). 

The  point  of  this  is  that  Equations  (1)  are  used  throughout  and  so  there  is  less  logic  to 
code  by  avoiding  a  zone  boundary  (especially  so  since  we  in  fact  use  <j>  and  rp  as  variables). 
Moreover,  because  an  implicit  solver  is  used,  information  is  spread  throughout  without 
any  lagging  of  zone  boundaries.  Now,  in  fact,  it  would  take  fewer  operations  per  iteration 
if  Equations  (1)  were  turned  off  in  the  conservation  law  zone.  But  the  code  would,  as  just 
stated  above,  be  more  complex  and  likely  require  more  iterations  between  zones  to  reach 
a  steady  state  if  this  were  done. 

What  we  have  done  in  this  zonal  method  is  to  use  a  single  “simple*  equation  set 
throughout.  In  “complex*  flow  zones,  a  r^re  complete  set  of  equations  is  used,  but  their 
effect  is  imposed  by  means  of  a  right  hand  side  forcing  function,  not  by  means  of  a  zonal 
boundary.  This,  I  believe,  makes  this  new  zonal  approach  much  more  versatile.  For 
example,  in  viscous  flow  we  should  be  able  to  evaluate  the  w  of  Equation  (lb)  from  a 
Navier-Stokes  or  boundary  layer  equation  zone.  Ideally  we  could  use  Equation  Set  (1) 
everywhere  and  use  more  complicated  equations  only  as  a  means  to  overwrite  w  and  a. 
Such  an  approach  is  currently  being  formulated. 

To  conclude,  my  feeling  is  that  zonal  codes  are  needed  and  have  their  place,  generally 
in  a  large  engineering  design  environment.  Small  engineering  teams  that  have  to  model  a 
variety  of  different  flow  fields  and  configurations  should  use  more  general  codes  that  are 
easier  to  set  up  for  a  given  problem.  In  designing  new  zonal  codes  it  is  important  that  we 
try  to  keep  the  number  of  equation  sets  and  zone  interface  boundaries  to  a  minimum  for 
simplicity.  The  type  of  zonal  approach  developed  in  Appendix  E  is  one  such  attempt,  and 
it  appears  to  be  very  promising. 


SECTION  m 
CONCLUSIONS 

Grid  generation  procedures  and  zonal  solution  methods  were  studied.  Both  elliptic  and 
hyperbolic  three  dimensional  grid  generation  procedures  were  developed.  The  hyperbolic 
grid  generator  is  especially  easy  to  use  and  it  is  very  efficient.  It  can  fail,  however,  whenever 
the  body  surface  is  discontinuous  or  the  used  specified  surface  grid  distribution  is  too 
irregular.  The  elliptic  solver  is  more  robust,  but  it  requires  much  more  user  input  and 
much  more  computer  time. 

FYom  our  experience  with  zonal  codes,  we  concluded  that  they  have  their  place,  partic¬ 
ularly  for  design  applications  where  maximum  computational  efficiency  is  requiied.  How¬ 
ever,  a  major  effort  must  be  made  to  reduce  the  complexity  of  zonal  codes.  An  approach 
which  interfaces  the  zones  through  forcing  functions  rather  than  boundary  conditions  was 
developed,  and  it  appears  to  offer  this  possibility. 
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ABSTRACT 

An  algorithm  for  generating  computational  grlda 
about  arbitrary  three-dimensional  bodies  la  developed. 

The  elliptic  partial  differential  equation  (FDE)  approach 
developed  by  Steger  and  Sorenaon  and  used  In  the  NASA  ' 
computer  program  CRAPE  Is  extended  from  two  to  three 
dimensions.  Forcing  functions  which  are  found  automati¬ 
cally  by  the  algorithm  give  the  user  the  ability  to  con¬ 
trol  mesh  cell  size  and  skewness  at  boundary  surfaces. 
This  algorithm,  as  Is  typical  of  PDE  grid  generators, 
gives  smooth  grid  lines  and  spacing  In  the  Interior  of 
the  grid.  The  method  is  applied  to  a  rectilinear  wind- 
tunnel  case  and  to  two  body  shapes  in  spherical 
coordinates. 

NOMENCLATURE 

a  coefficient  In  forcing  function  Influencing  decay 
rate  of  control  at  boundary 

b  coefficient  in  forcing  function  Influencing  decay 
rate  of  control  at  boundary 

c  coefficient  In  forcing  function  Influencing  decay 
rate  of  control  at  boundary 

J  Jacobian  of  transformation,  determinant  of  M 
M  matrix  of  transformation  metrics 
P  forcing  function  In  Poisson  equation 

P  factor  In  P.  giving  control  of  cell  size  and 
skewness  at  boundary 

Q  forcing  function  in  Poisson  equation 

Q  factor  In  Q,  giving  control  of  cell  size  and 
skewness  at  boundary 

R  forcing  function  in  Poisson  equation 

ft  factor  in  R,  giving  concrol  of  cell  size  and 
skewness  at  boundary 

r  column  vector  having  elements  x,y,z 

RHS  term  used  in  Eq .  (7)  In  solving  for  P.Q.R 

S  distance  along  line  of  Increasing  ;  In  real  domain 


T  superscript  indicating  transpose  of  a  matrix 

x  Independent  variable  in  real  domain,  Cartesian 
coordinate 

y  Independent  variable  In  real  domain,  Cartesian 
coordinate 

z  Independent  variable  In  real  domain,  Cartesian 
coordinate 

a  nonlinear  coefficient  in  Eq.  (2a) 

Y  signed  cofactor  of  the  matrix  M 

T  matrix  having  elements  y 

d  finite  difference 

C  Independent  variable  in  computational  domain 

n  Independent  variable  in  computational  domain 

8  angle  down  from  axis  toward  equator  In  spherical 
coordinates 

5  Independent  variable  In  computational  domain 

o  distance  from  origin  In  spherical  coordinates 
♦  angle  around  axis  In  spherical  coordinates 
ui  relaxation  parameter  for  polnt-SOR 

1  row  number  in  matrices  M  and  T 

j  column  number  In  macrlces  M  and  T 

INTRODUCTION 

The  ability  to  generate  grids  about  arbitrary  three- 
dimensional  aerodynamic  configurations  stands  today  as 
one  of  the  critical  pacing  items  In  computational  fluid 
dynamics  U).  Of  the  various  methods  for  generating 
grids,  the  elliptic  partial  differential  equation  tech¬ 
nique  (2-6) ,  with  Its  Inherent  smoothness,  has  proven  to 
be  one  of  the  most  automatic  and  general  approaches.  In 
two-dimensional  applications  this  has  been  especially 
true  when  the  elliptic  grid-generation  equations  are 
combined  with  an  algorithm  that  automatically  chooses 
inhomogeneous  terms  to  give  the  user  control  of  mesh 
cell  size  and  skewness  at  boundaries.  Such  an  algorithm, 
as  developed  in  Refs.  7  and  8,  has  resulted  in  the 
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widely  used  NASA  computer  program  GRAPE1  (8).  The  exten¬ 
sion  of  the  GRAPE  algorithm  to  three  dimensions  has  been 
long  overdue  and  Is  the  subject  of  this  paper. 

In  the  algorithm  developed  here,  three-dimensional 
elliptic  partial  differential  grid-generation  equations 
are  modified  co  give  Che  user  concrol  of  the  spacing  and 
skewness  of  mesh  lines  that  approach  the  boundary.  This 
Is  accomplished  by  adding  co  Che  grid  generation  equa¬ 
tions  a  sec  of  forcing  terms  which  are  automatically 
evaluated  by  imposing  differential  equation  constraints 
of  arc  length  control  and  surface  orthogonality  at  Che 
boundary.  The  user  has  only  to  Input  Che  desired  grid 
spacing  at  the  boundary  surface.  A  description  of  these 
grid-generation  equations  and  a  numerical  solution  pro¬ 
cedure  are  developed  In  Che  main  part  of  this  paper. 
Because  a  spherical  topology  is  one  of  the  most  efficient 
grid  systems  In  three  dimensions,  the  grid-generation 
algorithm  Is  further  modified  to  avoid  difficulties  in 
generating  grids  near  a  spherical  or  cylindrical  axis 
singularity.  These  details  and  grid  results  in  both 
rectangular  and  spherical  grid  topologies  are  presented 
in  the  remainder  of  the  paper. 

GRID-GENERATION  EQUATIONS 


The  elliptic  partial  differential  equation  grid- 
generation  approach  has  been  extensively  developed  else¬ 
where  (9).  Elements  of  this  approach  are  reviewed  below 
co  develop  source  cerms  chat  automatically  enforce 
interior  mesh  line  orthogonality  to  a  boundary  surface 
and  give  the  user  control  of  the  step  size  between  the 
boundary  and  Che  next  Interior  surface. 

Ic  is  required  that  the  mapping  between  physical 
space  x,y,z  and  computational  space  C.b.4  satisfy 
the  Poisson  equations 
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Given  proper  choice  of  the  source  terms  P,Q,R,  these 
equations  satisfy  the  maximum  principle  and  thus  ensure 
a  one-to-one  mapping.  Equation  (1)  is  conveniently 
solved  numerically  in  Che  uniform  computational  space, 
4,n,4.  The  equations,  so  transformed  (10-13) ,  are 


LU 


where 


a22rnn  +  a3ir44  +  2(ai*r4n  +  a”rU  +  a”rn4 


-J2(P?C  +  Q?n  +  R??)  (2a) 
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Ylj  Is  the  Ijth  signed  cofactor  of  the  matrix  M 

x.  > 


M 


'5  'n  '4 
z,  z  z 


(2d) 


and  the  Jacobian 


5  "n 

J  Is  the  determinant  of  M. 


In  the  present  application  Che  Inhomogeneous  terms 
P,Q,R  control  grid  spacing  and  skewness  for  mesh  cells 
adjacent  co  the  body  boundary  surface,  4  •  0.  The 
forcing  terms  are  chosen  to  be 

P(C.n.C)  -  P(S,n)e‘a;  (3a) 

Q(C.n,4)  -  Q«.n)e'b5  (3b) 

R(E,n,C)  -  R(C,n)e-c^  (3c) 

when  the  exponential  factors  cause  the  control  to  dissi¬ 
pate  or  relax  with  distance  from  the  boundary  4  «  0. 
Relaxing  the  control  with  distance  from  the  4  •  0 
boundary  is  necessary  so  as  not  to  overly  constrain  the 
grid  lines  with  respect  to  the  opposing  (4  -  tg^,) 
boundary.  The  positive  constants  a,b,c  Influence  the 
rate  of  decay  of  Che  boundary  control. 

In  most  elliptic  grid-generation  techniques,  the 
points  on  boundaries  are  user-specified.  Thus,  the  C,n 
distribution  of  grid  polncs  Is  specified  on  the  4  *  0 
surface.  Figure  1  shows  a  typical  mesh  cell  touching  a 


Fig.  1  Typical  mesh  cell  touching  4  •  0  surface 

4  ■  0  surface,  with  4  and  n  varying  over  that  surface. 
The  condition  for  orthogonality  of  the  grid  lines  inter¬ 
secting  the  4  -  0  boundary  surface  is  that  the  unit 
vectors  in  the  4  and  4  directions  and  in  the  n  and  4 
directions  be  mutually  normal.  These  two  conditions  can 
be  expressed  by  the  vector  dot  products 


rc  ■  r(  -  0  (4a) 

rn  •  r;  -  0  (4b) 


To  control  the  cell  size  on  the  t  ■  0  boundary,  the 
"height”  of  the  cells  must  be  regulated  to  some  pre¬ 
specified  value.  Letting  S  be  that  height,  the  dis¬ 
tance  along  a  line  of  increasing  ;,  we  wish  to  specify 
4S/44  at  the  4*0  boundary.  In  differential  form 
this  third  boundary-control  equation  can  be  expressed  as 


(4c) 


The  solution  procedure  described  below  requires  that 
Eq.  (4)  be  solved  for  the  derivatives  with  respect  co 
4  at  the  surface,  giving 


z 


4 


_ 3S/34 _ 

(Y13/y)3  +  Y2  j/Yj  J  +  Dl/J 


(5a) 


‘An  acronym  derived  from  "GRlds  about  Airfolla 
using  Poisson's  Equation." 


(5b) 


13 


1 


The  surface  control  aquations  expressed  by  althar 
Eqs.  (4)  or  (S)  ara  solved  simultaneously  with  tha 
Interior  grid-ganeration  equations  given  by  Eq.  (2). 
Because  x.y.z  are  specified  on  the  boundary,  the  sur¬ 
face  control  equations  supply  three  additional  relations 
that  can  be  used  to  determine  the  unknown  values  of 
P.Q.R.  An  Iterative  solution  process  Is  used  as 
described  below. 

SOLUTION  PROCEDURE 


An  Iterative  procedure  la  used  to  solve  the  grid- 
generation  equations.  Each  lteradon  is  In  two  distinct 
parts.  In  the  first  part  we  begin  with  the  x.y.z  from 
the  Initial  conditions  or  the  previous  Iteration,  and 
update  values  of  the  P.Q.R  terms.  The  second  part  of 
the  Iteration  step  updates  values  of  x.y.z  at  each 
point  In  the  field  using  the  new  values  for  P.Q.R. 

To  update  the  P.Q.R,  we  first  note  that  solving 
Eq.  (2)  produces  a  grid  for  an  appropriate  choice  of 
P.Q.R.  He  wish  to  Impose  constraints  on  grid  cells  at 
the  boundary  4  •  0,  and  thus  determine  P.Q.R.  Equa¬ 
tion  (4)  gives  the  constraints  and  Eq.  (5)  gives  the 
same  constraints  expressed  as  requirements  on  derivatives 
of  r  with  respect  to  4  at  the  boundary.  The  deriva¬ 
tives  In  Eq.  (5),  along  with  difference  approximations  ^ 
for  all  other  first  and  second  partial  derivatives  of  r 
with  respect  to  4,n,4,  are  substituted  Into  Eq.  (2), 
which  is  solved  for  P.Q.R. 

To  obtain  the  difference  approximations  for  first 
and  second  partial  derlvltives  of  r  with  respect  to  . 
$.0.4  on  the  4-0  boundary  surface,  we  proceed  as 
follows.  Because  x.y.z  are  specified  at  4-0.  first 
and  second  partial  derivatives  of  r  with  respect  to 
4  and  n  on  that  surface  can  be  found  by  differencing 
fixed  boundary  points.  Those  derivatives,  combined  with 
user  specification  of  ^S/34.  are  used  in  Eq.  (5)  to 
determine  derivatives  rr.  Second  partial  derivatives 
r^  and  rn(;  are  found  b\f  differencing  r.  with  respect 
to  4  and  n.  Thus,  the  only  derivatives  lacking  in 
Eq.  (2)  are  r_n.  These  are  found  by  differencing  the 
solution  for  r  at  and  near  the  surface  using  the  cur¬ 
rent  interior  grid  solution. 

Thus,  at  each  4-0  point,  Eq.  (2a)  are  three 
equations  which  can  be  solved  for  the  three  unknowns 
P.Q.R.  From  Eq.  (3)  It  can  be  seen  chat  at  4-0, 
P(4.n,4)  reduces  to  P(4.n).  and  similarly  for  Q  tnd  R. 
From  Eq.  (2)  RH?  is  defined  as 


52 


-J 


[Vr44  +  *”rnn  +  a>ir44 


+  2<“iir4n  +  aur44  +  °2  irn4)  1 

The  solution  then  is 


(6) 


M~l  RHS 


rT  rh?/j 


(7) 


where  T  Is  the  matrix  having  elements  Yij.  For  the 
whole  field  P.Q.R  are  then  found  by  multiplying  P,§,R 
by  the  appropriate  exponential  factors  as  in  Eq.  (3). 

The  second  part  of  each  iteration  step  is  to  use  the 
new  values  for  P.Q.R  in  Eq.  (2)  to  find  new  x.y.z 
everywhere  in  the  field.  In  this  research  effort  the 
Iterative  solution  procedure,  chosen  for  ease  of  coding, 
waa  polnt-SOR.  Thus ,  the  P.Q.R  terms  necessary  to 
cause  the  grid  to  have  the  desired  behavior  at  the  bound¬ 
ary  ara  found  automatically  along  with  the  x.y.z  in  the 


interior  as  this  two-part  Iteration  schema  proceeds  to 
convergence. 

SPHERICAL  GRIDS 

From  a  computational  point  of  view  a  spherical 
topology  Is  one  of  the  moat  efficient  grid  systems  for 
three-dimensional  bodies  because  a  spherical  grid  saves 
points  In  che  far  field.  For  example,  in  Navler-Stokes 
calculations  chat  use  highly  clustered  grids  naar  tha 
body  boundary  with  a  single  rectilinear  coordinate  sys¬ 
tem,  the  fine  grid  near  che  body  can  extend  Into  the  far 
field.  This  does  not  occur  with  spherical  coordinates. 

A  spherical  coordinate  system  does  Introduce  en  axis 
singularity,  but  experience  shows  chat  an  axis  singular¬ 
ity  can  be  readily  handled  in  flow  codes  that  use  general 
4.n,4  coordinates  and  solve  the  flow  equations  In 
conservative  form  (14) . 

For  grid  generation  using  Eq.  (2),  the  axis  singu¬ 
larity  of  the  spherical  coordinate  system  has  proven  to 
be  difficult.  Following  a  suggestion  of  J.  K.  Hodge 
(private  communication,  1977),  M.  Vinokur  and  J.  L. 

Sceger  (private  communication.  1978)  found  chat  one  wey 
to  avoid  the  axis  singularity  was  to  "convert"  Eq.  (1) 
into  a  pseudospherlcal  equation.  Instead  of  using  true 
spherical  Poisson  equations,  the  independent  variables 
x.y.z  are  simply  replaced  with  in  Eq.  (2)  for 

spherical  coordinates  or  by  p,$,z  for  cylindrical 
coordinates.  See  Fig.  2.  Thus,  Eq.  (1)  Is  replaced  by 

«pp  +  5«e  +  ‘  P(c'n>5)  (8*> 

npo  +  n68  +  %♦  ’  Qa,n,;)  (8b> 

So  +  ;99 +  S*  ‘  R(5>n>°  (8C> 


Fig.  2  Cut  away  sketch  showing  computational  variables 
and  spherical  coordinates 

Equation  (8)  is  not  in  the  form  of  a  true  Laplacian 
operator  in  spherical  coordinates.  It  is,  however,  an 
elliptic  equation  that  satisfies  che  maximum  principle 
and  it  generates  a  smoothly  varying  grid  just  as  Eq.  (1) 
does.  Ac  the  axis,  f  varies  monoconically  with  4  (see 
Fig.  2)  and  no  singular  behavior  is  encountered.  Trans¬ 
formation  of  Eq.  (8)  to  uniform  computational  space 
results  in  Eq.  (.’)  with  r  replaced  by 
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Solution  of  the  spherical  variable  fora  of  Eq.  (2) 
proceeds  such  aa  before.  Boundary-point  valuae  of  x,y,z 
are  apeclfied  along  with  an  initial  guaaa  of  the 
interior-point  values  of  x.y.z.  At  each  of  the  polnta, 
values  of  x.y.z  are  then  converted  in  the  usual  way  to 
spherical  values  of  0,8,1.  The  grid-generation  equa¬ 
tions  are  than  solved  by  relaxation,  and  once  a  solution 
is  obtained,  it  is  converted  back  to  x.y.z  values  at 
each  grid  point.  Equation  (8)  does  not  generate  the 
ease  solution  as  Eq.  (1),  but  the  solution  is  one  that 
can  be  Just  as  satisfactory.  Ic  must  be  cautioned  that 
if  1  varies  from  0  co  2ir  in  a  periodic  grid,  then 
one  must  difference  in  a  manner  which  accounts  for  a 
discontinuity  of  the  function,  but  not  of  the  derivative 
by  adding  or  subtracting  2it  where  appropriate. 

The  surface  clustering  and  orthogonality  relations 
aa  developed  previously  are  written  in  terns  of  x.y.z. 
Rather  chan  attempt  co  rework  these  relations  in  spheri¬ 
cal  variables,  we  simply  convert  p,6,$  variables  back 
co  x.y.z  variables  co  enforce  Eq.  (A).  That  is,  when¬ 
ever  P.Q.R  must  be  evaluated,  grid-point  values  of 
P,8,d  are  converted  back  co  x.y.z.  Values  of  P.Q.R 
are  found,  and  Chen  x.y.z  values  are  converted  back  to 
p.8,4.  Since  this  transformation  and  its  inverse  need 
be  done  only  for  Che  flrsc  few  (typically  3)  shells  of 
points  near  the  body  surface,  the  increase  in  computa¬ 
tion  is  not  significant. 


RESULTS 


The  method  has  been  coded  co  control  surface  orthog¬ 
onality  and  step-size  spacing  only  at  the  inner  (5  -  0) 
boundary.  The  same  3-D  computer  code  can  be  used  to 
generate  grids  either  for  rectangular  grid  topologies  • 
using  Eq.  (1)  in  terms  of  x.y.z  coordinates,  or  warped, 
spherical  grid  topologies  using  Eq.  (8)  for  p,8,$ 
coordinates. 


Rectilinear  Grids 


Rectilinear  three-dimensional  grids  with  automatic 
surface  step-size  concrol  have  been  generated.  A  rec¬ 
tangular  wind  tunnel  with  a  bump  on  the  floor  is  used  as 
a  test  case.  Since  this  case  fits  within  the  rectilinear 
topology,  Eq.  (1)  is  used  as  the  grid  generation  equa¬ 
tion.  Figure  3(a)  Illustrates  the  floor  (C  -  0)  boundary 
surface.  The  grid  in  this  case  has  40  points  in  the  E 
or  screamwlse  direction,  20  points  in  the  n  or  spanwlse 
direction,  and  23  points  in  the  q  or  vertical  direc¬ 
tion.  The  mesh  lines  Intersecting  the  floor  are  required 
co  do  so  normally,  and  a  spacing  normal  co  the  floor 
equal  to  0.0025  times  the  height  of  the  tunnel  is 
enforced.  That  spacing  is  approximately  l/20ch  of  what 
would  result  if  Che  points  were  equally  spaced  in  the 
vertical  direction. 

Figure  3(b)  shows  the  fifth  q  *  constant  coordi¬ 
nate  surface  above  the  floor.  The  control  of  point  spac¬ 
ing  normal  to  the  5  «  0  boundary  surface  is  shown  here 
by  the  uniform  proximity  of  this  surface  to  the  floor. 
Figure  3(c)  shows  Che  fifteenth  q  *  constant  surface, 
above  the  floor.  Here,  the  bump  has  begun  to  "flatten" 
since  the  grid  has  a  flat  upper-boundary  plane.  An 
example  of  a  different  family  of  coordinate  surfaces,  a 
E  ■  constant  surface  is  shown  in  Fig.  3(d).  ft  rides 
the  crest  of  the  bump  from  side  to  side  in  the  tunnel. 
Figure  3(e)  shows  an  example  of  the  third  family  of 
coordinate  surfaces,  an  n  -  constant  surface.  It  is 
the  fifth  such  surface,  counting  from  the  near  side  of 
the  tunnel.  The  viewpoint  here  la  normal  to  that  sur¬ 
face.  Figure  3(f)  shows  a  closeup  of  a  region  of  the 
surface  shown  on  the  previous  Fig.  3fe>,  It  Is  the 
region  at  the  right-hand  side  of  the  base  of  the  bump. 
Note  that  the  spacing  normal  to  the  bottom  boundary  is 


£11111 


Fig.  3  Grid  for  rectangular  wind  tunnel  with  bump  on 

floor 

fa)  Inner  (q  -  0)  boundary  surface  (rhe  floor), 
showing  bump 

(b)  Fifth  q  *  constant  coordinate  surface  above 
the  floor 

(c)  Fifteenth  q  -  constant  coordinate  surface 
above  the  floor 
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Fl«.  3  Concluded 

(d)  C  •  constant  coordinate  surface  on  top  of 
bump 

(e)  Fifth  o  •  constant  surface  passing  over 
top  of  the  buap.  orthogonal  view 

(f)  Closeup  of  region  near  base  of  bump 


constant  and  very  email.  That  spacing  would  tend  to  be 
much  larger  and  nonunifora  In  this  concave  region  If  Che 
grid  had  been  generated  without  the  forcing  teraa. 

Warped  Spherical  Grids 

Grids  for  fuselage  shapes  and  wing  shapes  have  been 
generated  using  the  spherical  coordinate  systsa  and 
solving  Eq.  (8).  As  of  this  writing  soaa  convergence 
pr  <bleas  are  encountered  with  fins  grids  using  spherical 
topology.  This  Is  currently  under  Investigation.  The 
results  presented  here  for  spherical  topology  are  Halted 
to  moderately  coarse  grids. 

A  grid  having  20  points  in  the  clrcuaferentlal  (O 
direction.  19  pointa  In  the  pole-to-pole  (n)  direction, 
and  30  points  In  the  outward  (()  direction  has  been 
generated  about  a  10  by  2  by  1  ellipsoid.  Figure  4(a) 
Illustrates  the  specified  5  and  n  distributions  about 
tha  top  half  of  this  ellipsoidal  body.  Figure  4(b) 
shows  the  n  ■  constant  coordinate  surface  about  the 
"equator,"  while  Fig.  4(c)  shows  s  closeup  of  the 
10  Innermost  lines  in  the  equatorial  plane.  Note  that 
the  spacing  normal  to  the  body  is  controlled.  Fig¬ 
ure  4(d)  shows  a  different  n  -  constant  coordinate 
surface,  this  one  nearer  to  the  pole  than  the  equator. 

It  Is  not  s  plane,  but  rather  more  In  the  shape  of  a 
bowl,  open  toward  the  lower  right.  An  example  of  a  dif¬ 
ferent  family  of  coordinate  surfaces,  a  £  •  constant 
surface,  Is  shown  In  Fig.  4(e),  extending  froa  pole  to 
pole.  Figure  4(f)  Is  a  closeup  of  the  region  of 
Fig.  4(e)  near  the  left-hand  pole  and  near  the  body.  The 
results  Indicate  that  the  surface  orthogonality  control 
and  the  clustering  control  are  working  well,  and  that  the 
solution  about  the  axis  is  well  behaved.  Examples  of  the 
third  kind  of  coordinate  surface,  4  •  constant,  are 
visually  Indistinguishable  from  the  body,  Fig.  4(a). 

On  an  axis,  e  Is  either  0  or  *  and  values  of  P 
and  o  are  obtained  by  parabolic  extrapolation  to  the 
axis  from  the  two  nearest  points  In  the  n  direction. 

The  three  conditions  used  to  determine  the  parabola  are 
that  the  mesh  line  must  pass  through  the  two  nearest 
points  111  the  r  direction,  and  that  It  must  have  zero 
slope  relsclve  to  the  axis  ss  It  crosses  the  axis. 

Values  of  for  each  point  on  the  axis  could  be  found 

by  extrapolating  along  any  of  the  lines  encountered  as 
4  varies  trom  0  to  (ma*.  and  in  fact  the  p  used  is 
an  average  of  all  those  values. 

Some  views  ot  another  spherical  grid,  this  one 
having  30  by  19  by  30  points,  are  shown  In  Fig.  5.  This 
grid  Is  about  a  wing  having  the  same  5  to  1  ratio  of  span 
to  midchord  elliptical  planform  as  the  ellipsoid  above, 
and  a  1 9Z  thick  airfoil  section.  The  only  adaptation 
necessary  here  is  that  the  P.Q.R  terms  at  the  trailing 
edge  sre  replaced  by  the  average  of  their  values  Imme¬ 
diately  above  and  below  the  edge,  l.e.,  they  sre  "aver¬ 
aged"  across  the  edge  In  the  E.  direction.  The  airfoil 
section  and  the  Innermost  four  grid  lines  are  shown  in 
Fig.  3(a).  The  norms)  spacing  and  the  angularity  of  the 
lines  intersecting  t lie  oidv  ate  nicely  controlled.  The 
sharp  trailing  edge  is  treated  succesatully,  as  shown  in 
Fig.  3(b).  Figure  i’ci  shows  the  far-fleld  behind  the 
trailing  edo.  ,  and  illustrates  the  ability  of  the  PDE 
method  to  generate  a  smooth  grid  over  a  sharp  corner. 

In  all  of  the  above  cases,  values  for  the  relaxation 
parameter  .  In  the  poinc-SOR  varied  from  0.8  to  1.6. 
Values  ot  the  a.o.c  parameters  In  Eq.  (3)  were  approxi¬ 
mately  0.3.  The  lumper  of  Iterations  necessary  varied 
from  25  tc  150. 


If) 


* 


FI,.  4  Grid  about  ellipsoid  ualng  spherical  coordinates  Fig.  4  Concluded 


(a) 

(b) 


(c) 


Top  half  of  ellipsoidal  body  surface 
Tenth  n  •  constant  coordinate  surface, 
l.e.,  equatorial  plane 
Cloaeup  of  lnnernoat  10  clrcuaferentlal 
lines  In  equatorial  plane 


Cd) 


(e) 


(f) 


17 


Fifth  bowl  (n  -  constant  coordinate  sur¬ 
face),  counting  toward  equator  frost  pole 
5  *  constant  coordinate  surface  extending 
frost  pole  to  pole 

Closeup  of  C  *  constant  coordinate  sur¬ 
face,  near  pole  and  near  body 


■i  1  ■  tii  in.  ■  iiin  s’l-eks aAteftew- — 


CONCLUSION 


Grid  abouc  wing  having  airfoil  section  and 
elliptical  planform 

(a)  Closeup  of  innermost  five  circumferential 
lines  in  equatorial  plane 

(b)  Closeup  of  trailing  edge  in  equatorial  plane 

(c)  Far-field  behind  trailing  edge 


An  elliptic  partial  differential  grid-generation 
technique  which  gives  the  user  ability  to  control  mesh 
cell  size  and  skewness  at  a  boundary  has  been  general¬ 
ized  from  two  to  three  dimensions.  The  method  discussed 
in  this  paper  has  been  applied  successfully  to  a  variety 
of  topologies  and  teat  cases.  Future  plana  Include  sub¬ 
stituting  a  faster  solution  method  such  as  ADI,  investi¬ 
gating  the  application  of  the  method  to  a  wider  collec¬ 
tion  of  topologies,  and  writing  a  transportable ,  user- 
oriented  three-dimensional  code,  such  as  GRAPE  is  In 
two  dimensions. 
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APPENDIX  B 

GENERATION  OF  THREE  DIMENSIONAL  BODY  FITTED  COORDINATES 
USING  HYPERBOLIC  PARTIAL  DIFFERENTIAL  EQUATIONS 


INTRODUCTION 


Body  conforming  curvilinear  grids  are  often  used  in  finite  difference  flow  field  simula¬ 
tions.  One  reason  for  this  is  that  the  application  of  boundary  conditions  can  be  simplified 
in  finite  difference  calculations  because  grid  lines  coincide  with  boundary  lines.  This  is 
especially  important  in  high  Reynolds  number  viscous  flow  simulation  in  which  high  flow 
gradients  near  the  body  surface  must  be  resolved. 

The  task  of  generating  a  satisfactory  body  conforming  coordinate  system  is  not  easy. 
The  grids  must  not  be  too  distorted,  they  should  have  smooth  variation,  and  they  should 
be  clustered  to  flow  field  action  regions  -  typically  near  boundary  surfaces.  Moreover,  the 
grids  should  be  generated  in  an  automatic  manner  that  requires  a  minimum  of  user  input. 

One  approach  for  generating  body  conforming  grids  with  minimum  user  input  has 
been  to  solve  a  set  of  partial  differential  equations.  In  this  technique  level  lines  of  £(x,  y,  z), 
r)(x,y,  z),  and  f(z,  j /,  z)  that  have  monotone  variation  are  sought  as  a  solution  of  a  Bet  of 
partial  differential  equations.  Generally  values  of  £,  f?  and  f  are  user  specified  on  the 
boundary  surface  and  constraints  expressed  as  differential  equations  are  used  to  develop 
the  grid  away  from  the  boundaries.  The  most  popular  such  approach  requires  the  solution 
of  a  set  of  elliptic  equations  that  satisfy  the  maximum  principle,  however,  hyperbolic 
and  parabolic  governing  equations  have  been  used  as  well,  at  least  in  two  dimensional 
applications. 

In  this  appendix  one  way  of  extending  the  hyperbolic  grid  generation  method  of  Steger 
and  Chaussee  to  three  dimensions  is  developed.  In  two  dimensions  the  two  differential 
constraints 


6*7.  +  =  0 

(la) 

=  (AVT1 

(lb) 
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Iff 
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or  in  £,  rj  computational  space 


x(xn  +  VtVn  =  0  (2o) 

x(yn  ~  Wi  =  AV  (26) 

have  been  solved  by  marching  in  rj  from  an  initial  data  plane  rj(x,y)  =  constant.  The  first 
equation  is  a  constraint  of  orthogonality.  The  second  equation  controls  the  mesh  spacing 
with  the  user  specifying  the  mesh  control  volume  AK  (actually  area  in  two  dimensions). 
A  linearized  version  of  equations  (2)  is  readily  shown  to  be  hyperbolic  and  suitable  for 
marching  in  rj.  Equations  (2)  are  solved  in  computational  space  to  give  the  x,y  location 
of  the  £  =  constant  and  17  =  constant  grid  lines. 

The  two  partial  differential  equation  ,  .  pressed  as  either  Equations  (1)  or  Equations 
(2),  have  been  referred  to  as  a  mesh  cell  volume  procedure  for  grid  generation.  In  the  next 
section  a  three  dimensional  procedure  is  developed. 


THREE  DIMENSIONAL  GRID  GENERATION  EQUATIONS 


A  body  fitted  exterior  grid  about  an  arbitrary  closed  boundary  surface  is  desired. 
Only  a  simple  topology  such  as  that  illustrated  in  Figure  (1)  will  be  considered  here.  The 
body  surface  is  chosen  to  coincide  with  c(x ,  y,z)  —  0  and  the  surface  grid  line  distributions 
of  (  =  constant  and  t)  ~  constant  are  user  specified.  The  outer  boundary  f(x,  y,z)  =  fm0, 
is  not  specified,  it  is  only  required  to  be  sufficiently  far  removed  from  the  inner  boundary. 
Using  f  as  the  marching  direction,  partial  differential  equations  are  sought  which  produce 
planes  of  constant  and  f  to  form  a  nonsingular  mesh  system. 

An  extension  of  the  mesh  cell  volume  procedure  to  three  dimensions  is  proposed.  In 
three  dimensions,  however,  there  are  three  orthogonality  relations  and  one  cell  volume 
constraint.  At  any  point  four  equations  n—  available  to  predict  the  three  unknowns  x,y 
and  z  so  one  equation  must  be  discarded.  Because  f  is  the  marching  direction  it  is  natural 
to  use  only  the  two  orthogonality  relations  that  involve  f,  this  leads  to  the  governing 
equations 


*{*{  +  ViV,  +  z<2(  o 


xnx,  4  y fjV<  ^  "  t) 


X(y„Z(  +  ZfS>(Zn  +  x*y(z,,  I,,y(z<  ~  x<y,,z(  =  AV 


or  with  f  defined  as  (x,  y,  z) 


u  r'« =■-  °-  fn  r'-  -  i  j  •-  •-  ^ 


(за) 

(зб) 

(3c) 


The  first  two  equations  represent  orthogonality  relations  between  {  and  ;  and  between  17 
and  f,  whiie  the  last  equation  is  the  volume  or  finite  Jacobian  constraint. 

Equations  (3)  comprise  a  system  of  nonlinear  partial  differential  equations  in  which 
x, y  and  z  are  specified  as  initial  data  at  —  0.  As  developed  below,  linearization  and 


analysis  of  Equations  (3)  about  a  nearby  known  state  reveals  that  the  system  is  hyperbolic 
with  f  as  the  marching  direction. 

Let  x°,  y°,  2°  represent  a  nearby  known  state  so  that 


*  =  i°  +  (x  —  i°)  =  x°  +  x 

y  =  y°  +  il  (4) 

Z  =  Z°  +  2 


where  x,y  and  z  are  small.  Substitution  of  these  expressions  into  Equations  (3)  and 
elimination  of  products  of  tilde  terms  results  in  the  locally  linearized  system 


-4o(r  —  r<j){  +  B0(r  —  r0)n  -+•  C0(r  —  f0)(  —  f 
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Let  R  =  r  -  r0,  then  (5)  is  rewritten  as 


A0R(  +  BqRt,  4-  C0R(  =  /  (7) 

Now  Cgl  exists  unless  (AV)_1  — »  oo,  which  we  will  not  impose,  so  (7)  can  be  rewritten  as 

CglAoR(  +  Ci'Boh  +  R(  =  C^J  (8) 

Although  the  verification  is  nontrivial,  Cq1A0  and  Cq1B0  are  found  to  be  symmetric 
matrices  (this  was  carried  out  by  Dennis  Jesperson  of  the  NASA  Ames  Research  Center, 
who  used  MACSYMA).  The  linearised  system  Equation  (8)  is  therefore  hyperbolic  and 
can  be  marched  with  f  serving  as  the  “lime -like”  direction. 

It  can  be  pointed  out  that  an  analysis  was  attempted  for  the  three  orthogonality 
relations  alone.  These  equations,  however,  are  readily  shown  to  be  improperly  posed  for 
marching  with  initial  data  in  Indeed,  as  best  as  we  can  discern,  the  three  relations 
do  not  lend  themselves  to  unique  solutions  regardless  of  the  type  of  boundary  conditions 
specified. 
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SOLUTION  PROCEDURE 


The  nonlinear  system  of  grid  generation  equations  given  by  Equations  (3)  are  solved 
with  a  noniterative  implicit  finite  difference  scheme.  An  unconditionally  stable  implicit 
scheme  is  chosen  so  the  marching  step  size  in  (  can  be  arbitrarily  selected  based  only  on 
considerations  of  accurately  generating  the  grid.  Iterative  solution  of  the  nonlinear  grid 
generation  equations  is  avoided  by  expanding  the  equations  about  the  previous  marching 
step.  As  a  consequence  Equation  (7)  is  solved  with  the  nearby  known  state  0  taken  from 
the  previous  f  step. 

a)  Numerical  Method 

Let  =  Arj  —  Af  =  1  such  that  £  =  j  —  l,tj  =  k  —  1  and  (  =  I  -  1.  Central  spatial 
differencing  of  Equations  (5)  in  £  and  r?  wisa  first  order  backward  implicit  differencing  in 
£  leads  to 

AiS((rt+1  -  ft)  +  BiSn(fi+l  -  rj)  +  C,  y?  n+i  =  9i+i  (9) 

where 


Vf?/+i  —  *!+ 1  —  ri 

Note  that  Q6tri  was  subtracted  from  /j+i  to  produce  gj+1  in  theh  above.  Throughout  only 
those  indices  that  change  are  indicated,  thus  r<+1  =>  =>  ry+li4,»  etc. 

Multiplying  through  by  Cf1  gives 

CrlA,S((f,+l  -  rj)  +  CrlB,6n(r,+l  -  r,)  +  J(fl+1  -  r()  =  C,-lg,+1  (10) 
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where  I  is  the  identifying  matrix.  To  reduce  the  inversion  cost  the  difference  equations  are 
approximately  factored  as 


(/  +  CrlA,6()(I  +  C,~l BiSn)(fi+i  -  f,)  =  Crlg,+l  (11) 

so  that  r/+i  is  obtained  by  solving  sequences  of  one-dimensional-like  block  tridiagonal 
systems 

(/  4-  C,  lAiS()gl+1  =  g/+l  (12a) 

(/  +  C^B^)  vf  f/+i  =  $+l  (126) 


r/+l  =  r,  +  Vfiv+,  (12c) 

Although  not  shown,  numerical  dissipation  terms  are  added  in  f  and  r)  directions. 

The  coefficient  matrices  At,  Bi  and  C/  contain  £  and  rj  derivatives  which  are  formed 
using  central  differences.  These  matrices  also  contain  derivatives  for  x(,y(  and  z(  which 
are  obtained  from  Equations  (3)  in  terms  of  £  and  T)  derivatives.  That  is,  Equations  (3) 
are  linear  in  the  unknowns  zf,yf  and  zf.  They  are  easily  solved  for  as 


y( 

Z' 


AV 

(DetC) 


/  x(z„  -  ynz(  \ 
x nz(  ~ 

\  Wn  ~  ) 


(13) 


with 


Det(C)  =  ( y(zn  -  ynz()2  +  {xnz(  -  x (zn)2  +  ( x(y„  -  xny()2 

Note  that  AV  /  +  y*  +  z*)  =  Det(C)  so  that  Det(C)  will  be  zero  if  and  only  if  the 

user  specified  AV  =  0.  Hence,  C~l  will  exist. 
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b)  Cell  Volume  Specification 


The  user  has  control  of  the  grid  by  means  of  the  initial  surface  point  distribution 
and  by  specification  of  the  cell  volumes,  Through  the  cell  volumes  the  extent  and 

clustering  of  the  grid  can  be  essentially  controlled.  Because  the  cell  volume  at  each  point 
must  be  specified,  it  is  clear  that  the  user  must  devise  some  kind  of  method  for  determining 
volumes.  There  are  many  possibilities,  here  one  such  approach  is  illustrated. 

Suppose  we  had  a  sphere  to  grid.  A  reasonable  grid  might  have  uniform  angle  spacing 
and  have  a  radial  grid  distribution  that  is  exponential.  For  this  special  geometry  and 
grid  we  can  analytically  determine  the  control  volumes  by  a  simple  formula.  Now  take 
the  problem  at  hand,  perhaps  an  aircraft  fuselage,  which  we  want  to  mesh  as  a  warped 
spherical-like  grid.  We  can  find  a  sphere  fbat  has  the  same  surface  area  as  our  fuselage 
and  use  the  grid  cell  volumes  of  the  sphere  to  specify  the  cell  volumes  of  the  fuselage  grid. 
However,  the  fuselage  will  not  have  the  same  kind  of  surface  area  distribution  as  a  sphere 
with  equal  angle  distribution.  So  here  we  need  an  adjustment,  something  like 


A Vjxi  -  {AVt  kj),phtre 


trt 


6 


(14) 


where  6  — ►  1  for  small  l  and  6  — »  0  for  large  l.  That  is,  the  volumes  would  be  adjusted 
initially  to  the  local  boundary  surface  increments.  But  as  we  march  out  the  uniform 
spherical  volumes  would  gradually  be  specified.  Such  an  approach  has  been  used,  and, 
as  a  result,  the  far  field  portion  of  the  grid  tends  to  be  uniformly  spherical.  The  results 
shown  in  the  next  section  illustrate  this  behavior. 
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RESULTS 


A  series  of  plots  are  shown  in  Figures  (1)  -  (4)  to  indicate  the  generation  of  grids 
about  ellipsoids  and  cambered  ellipsoid  test  geometries.  Figure  (1)  sets  the  grid  notation 
while  Figures  (2)  -  (4)  show  some  typical  results.  The  grids  shown  in  Figure  (2)  are  for 
an  ellipsoid  which  has  major  to  minor  axis  ratios  of  4  to  1  and  2  to  1.  The  views  2a  to  2c 
show  the  user  specified  surface  point  distributions  from  various  observer  positions.  Similar 
views  are  shown  in  Figures  (2d)  to  (2i)  after  marching  3  and  19  steps  in  the  f  direction. 
The  views  shown  in  Figure  (3)  show  a  grid  generated  about  a  cambered  ellipsoid.  Finally 
the  views  shown  in  Figure  (4)  show  a  grid  about  a  wing-like  ellipsoid  ratioed  6:1:  1/6. 
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Figure  2a.  Surface  Distribution  on 
View  at  17  =  0. 


Ellipsoid,  (  =  0. 
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Figure  2c.  Surface  Distribution  on  Ellipsoid 
View  from  £  =  0  “North  Pole.” 
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Figure  2f.  Grid  at  f  =  3Af. 

View  at  £  =  0. 
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Figure  2g.  Grid  at  f  =  19Af. 
View  at  t}  =  0. 
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Figure  3a.  Cambered  Ellipsoid,  f 


Figure  3b.  Grid  About  Cambered  Ellipsoid 
At  Surface  f  =  19Af. 


Figure  3c.  Grid  About  Cambered  Ellipsoid. 

Side  View,  Upper  Right  Quarter. 


Figure  3f'  v «;'id  AK.ut  Carubeici  L.upi  .jit. 

S>.Jc  V:-?w,  Upper  Left  Quarter. 
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Figure  4b.  Surface  Distribution  of  6  :  1  :  1/6  Ellipsoid 
Span  =  6,  Chord  =  1,  Thickness  =1/6 


Figure  4c.  North  Pole  View  of  6  :  1  :  1/6 
Surface  Distribution 
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Figure  4d.  Surface  Distribution,  Partial  Planform 
View  of  6  :  1  :  1/6 
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Figure  4g.  North  Pole  View  -  25th  Grid  Plane 
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ABSTRACT 


A  mesh  system  composed  of  multiple  overset  body- 
conforming  grids  is  described  for  adapting  finite- 
difference  procedures  to  complex  aircraft  configura¬ 
tions.  In  this  so-called  "chimera  mesh,"  a  major  grid 
is  generated  about  a  main  component  of  the  configuration 
and  overset  minor  grids  are  used  to  resolve  all  other 
features.  Methods  for  connecting  overset  multiple 
grids  and  modifications  of  flow-simulation  algorithms 
are  discussed.  Computational  tests  in  two  dimensions 
indicate  that  the  use  of  multiple  overset  grids  can 
simplify  the  task  of  grid  generation  without  an  adverse 
effect  on  flow-field  algorithms  and  computer  code 
complexity. 

INTRODUCTION 

Finite-difference  simulations  of  flow  about 
realistic  aircraft  configurations  are  likely  to  use 
more  than  one  grid  system  and  more  than  one  governing 
equation  set.  A  composite  computer  code  for  such  simu¬ 
lations  is  thus  chimera-like  —  a  chimera  being  a  mytho¬ 
logical  creature  with  the  head  of  a  lion,  the  body  of 
a  goat,  and  the  tall  of  a  snake. 

A  finite-difference  computer  code  that  uses  mul-  . 
tiplc  equation  sets  (e.g.,  zones  of  Navier-Stokes  and 
velocity  potential)  and  multiple  grids  will  be  consid¬ 
erably  more  complicated  than  a  code  that  uses  a  single 
flow-field  equation  set  and  a  single  grid  mapping. 
Nevertheless,  such  a  chimera  ode  will  lead  to  increased 
computational  efficiency  in  flow  simulation  of  complex 
aircraft  configurations.  As  a  result,  the  computational 
aerodynamic  1st  will  likely  spend  considerably  mure  time 
in  the  future  developing  ways  to  interface  grids,  gov¬ 
erning  equations,  and  data  bases. 

The  purpose  of  this  paper  is  to  explore  a  chimera- 
type  mesh  scheme  in  which  a  major  grid  is  generated 
about  a  main  body  clement  such  as  a  wlng-fusclagc,  and 
minor  grids  are  overset  on  the  major  grid  so  as  to 
resolve  secondary  features  of  the  configuration  such  as 
■tores,  nacelles,  and  canards.  In  general,  the  minor 
grids  are  overset  on  top  of  the  major  grid  without 


requiring  mesh  boundaries  to  join  in  any  special  way. 
With  the  use  of  such  an  overset  mesh  system  the  task 
of  grid  generation  is  simplified  because  individual 
grids  can  be  generated  almost  independently  and  then 
superimposed  to  form  the  overall  mesh  system.  In  this 
way,  one  can  build  up  a  grid  system  for  treating  com¬ 
plex  configurations  and  avoid  severe  mesh  distortions. 

Overset  mesh  systems  have  been  used  previously. 
The  explicit  finite-difference  code  of  Magnus  and 
Yoshihara  (1)  for  solving  the  transonic  flow  about  air 
foils  is  an  early  example  of  overset  grids  used  to 
achieve  computational  efficiency.  More  recently, 
Berger  and  Oliger  (_2)  have  developed  a  numerical  pro¬ 
cedure  which  automatically  inserts  overset  grids  to 
resolve  gradient  regions  of  two-dimensional  convection 
problems.  As  to  applications  to  complex  geometries, 
Starlus  (3)  has  used  such  an  approach  for  the  shallow- 
water  equations,  and,  in  work  begun  at  Ames  Research 
Center,  Atta  (4)  has  coded  overset  grids  for  the 
transonic  potential  equation.  Atta  has  also  extended 
the  concept  to  three  dimensions  (.5) .  We  also  have 
long  argued  the  potential  advantages  for  such  a  mesh 
system  (6,7) • 

In  the  research  described  here  we  have  undertaken 
an  independent  development.  This  is  in  order  to  adapt 
our  grid  system  to  an  existing  class  of  implicit 
finite-difference  algorithms  for  solving  the  unsteady 
potential,  Euler  and  thin-layer  Navier-Stokes,  and 
parabolized  Navier-Stokes  equations.  These  codes  all 
use  similar  coordinate  transformations,  and  zonal-flow 
equation  versions  are  under  way.  Ultimately  we  intend 
to  explore  the  complexities  of  a  true  chimera  code  by 
combining  the  various  equation  sets  with  the  overset 
grid  scheme. 

In  this  paper  we  develop  some  of  the  philosophy 
and  details  of  tin  overset  mesh  system.  The  methodol¬ 
ogy  here  Is  restricted  Lo  two  dimensions,  and  the  pres 
ent  solutions,  which  utilize  the  stream  function,  are 
intended  for  Rrid  evaluation.  Overset  grid  calcula¬ 
tions  using  the  Euler  and  thin-layer  Navier-Stokes 
equations  will  he  presented  elsewhere  in  a  subsequent 
publication  by  Benek,  Stcgcr,  and  Dougherty. 


APPROACH  AND  MOTIVATION 


We  wish  to  build  a  mesh  system  in  which  there  ie  a 
major  grid  generated  about  a  dominant  body.  Minor  grid 
systems  are  generated  about  remaining  portions  of  the 
body  boundary  or  are  used  to  better  resolve  portions  of 
the  major  grid.  At  this  point  we  do  not  allow  the 
minor  grids  to  communicate  with  each  other  except 
through  the  major  grid.  The  sketches  shown  in  Figs.  1 
to  4  are  presented  with  the  discussion  below  to  clarify 
these  ideas. 

Figure  1  illustrates  an  airfoil-flap  combination. 

A  major  grid  is  generated  about  the  main  airfoil;  a 
minor  grid  is  used  to  resolve  the  flap.  The  two  grids 
ire  not  joined  at  a  common  boundary  and  so  are  con¬ 
nected  as  follows:  data  from  the  major  grid  are  inter¬ 
polated  to  supply  outer  boundary  conditions  to  the 
minor  grid.  Thus,  the  flap  grid  receives  input  froei 
the  main  airfoil.  Within  some  curve  of  the  minor  grid 
that  circumscribes  the  flap,  points  of  the  major  grid 
will  be  blanked  out.  The  major  grid  points  forming  a 
perimeter  to  the  blanked  region  will  also  be  flagged. 
Flow  variables  at  these  perimeter  points  are  supplied 
by  interpolating  the  minor  grid  solution.  Thus  the 
effect  of  the  flap  is  Imposed  on  the  main  airfoil. 

The  grids  shown  '.a  Fig.  2  illustrate  another  appli¬ 
cation  of  overset  grids.  Here  a  main  rectangular-like 
grid  fits  well  to  a  cascade  blade  element  everywhere 
except  at  the  blunt  nose  region.  A  minor  mesh, 
wrapped  about  the  nose,  is  inserted  to  resolve  the  flow 
suction  peak.  The  edge  boundaries  of  the  minor  grid 
are  Interpolated  from  the  major  grid.  Points  in  the 
major  grid  in  the  vicinity  of  the  nose  are  blanked,  and 
the  flow  values  are  found  by  Interpolating  data  from 
the  fine  minor  grid  wrapped  around  the  nose.  Thia 
choice  of  overset  grids  can  be  preferable  to  the  use  of 
a  C-grid  throughout  because  a  very  skewed  grid  results 
for  blades  with  high  solidity,  as  indicated  in  Fig.  3. 


There  are  also  disadvantages  to  the  overset  mesh 
system.  Interpolation  points  and  blanked  points  must 
be  located  and  labeled  for  special  treatment.  The  solu¬ 
tion  algorithm  also  becomes  more  complex  relative  to 
using  a  single  grid  mapping.  Finally,  Interpolation 
can  cause  Inaccuracies  such  as  local  loss  of  conserva¬ 
tive  form. 


ALGORITHM  CHANGES 

To  use  overset  grids  the  finite-difference  algo¬ 
rithms  need  only  be  altered  In  three  esaentlal  ways. 
First,  the  data  base  must  be  structured  for  a  number  of 
grids,  each  of  which  can  have  a  different  dimension. 

This  is  easily  accomplished  In  a  variety  of  ways  depend¬ 
ing  on  the  computer.  For  example,  on  a  single 
instruction-stream  machine,  use  of  a  single  array  can 
be  practical.  In  this  case  the  index  j,k  of  the  nth 
grid  can  be  locaced  as  a  single  array-point  1.  The 
FORTRAN  variable  Q(l)  Is  located  by  the  index 

I  -  J  +  (K  -  1)  *  JMAX(M)  +  ISUM(M) 

where  JMAX(M)  and  K2iAX(M)  are  maximum  j  and  k  values 
of  the  Mth  grid  and  ISUM(M)  is  the  number  of  grid 
points  in  all  grids  before  the  Mth  grid.  Using  such 
a  single  Index  leads  to  a  dependency  which  does  not 
vectorize  on  some  machines.  In  which  case  other  con¬ 
structions  should  be  used. 


Second,  the  algorithm  must  skip  blanked  points. 

In  time-like  marching  algorithms  this  is  easy  to  Imple¬ 
ment  by  simply  multiplying  by  a  0  or  1  flag  array 
scored  at  each  point.  Such  blanking  is  illustrated 
below  for  the  centrally  differenced  Laplaclan  In  which 
an  Implicit  approximately  factored  delta-form  algorithm 
is  used  as  part  of  an  iterative  solution  process.  That 

is. 


Sxx*J,k  +  V".V 


(1) 


Figure  4  illustrates  multiple  overset  grids.  In 
this  case  the  minor  grid  resolving  the  airfoil  leading- 
edge  flap  (slat)  is  shown  intersecting  the  main  airfoil. 
Such  grid  points  must  be  turned  off  and  a  fringe  bound¬ 
ary  must  be  defined  about  these  points  with  data  sup¬ 
plied  by  Interpolating  the  major  grid. 


where 

4x**j  -  (Vi  ■  ui +  v»)/(Ax)1 
Vk  ’  (V>  "  2*k  +  *k-,)/(Ay>1 


(2) 


A  system  of  overset  grids  has  many  advantages.  It 
can  be  used  to  treat  complex  geometries,  resolve  large 
flow-field  gradients,  and  eliminate  grid  distortion. 
Overset  grids  ate  easier  to  construct  than  grids  gen¬ 
erated  by  patching  meshes  together  at  common  boundaries 
becauae  each  grid  is  somewhat  independently  generated. 
This  Independence  ensures  that  each  grid  will  maintain 
a  well-ordered  sec  of  points  that  are  successively 
indexed.  Well-otderednesa  is  an  important  property  of 
a  mesh  because  when  finite-difference  approxlmatione 
are  applied  to  it,  a  set  of  structured  algebraic  equa¬ 
tions  results.  The  computational  efficiencies  gained 
by  using,  for  example,  spectral  methods,  approximate 
factorization  or  alternating  directions  techniques,  and 
vectorized  computer  programming  require  structured 
matrices.  Another  advantage  of  overset  grids  is  a 
natural  occurrence  of  overlapping  grid  points  between 
Interpolated  boundaries.  If  the  grids  are  sufficiently 
well  overlapped,  the  stability  of  an  implicit  algorithm 
should  not  be  adversely  affected,  even  though  the 
interpolated  boundary -values  are  updated  in  an  explicit 
mode.  For  a  similar  reason,  iterative  convergence  can 
be  improved  by  overlapping  mesh  boundaries  (8). 


The  factored  difference  equations  can  be  represented  as 
(I  -  -  w6yy)(*"+1  -  *n)  -  0.(4^  +  Syy)*"  (3) 

or  la  algorltha  font  as 
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One  can  flag  off  points  by  simply 
“  "  =uJ.k  wJ,k  ’ 
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Alternately  one  can  Initially  proceed  as  If  none 
of  the  points  la  special.  Then,  just  before  solving 
the  first  trldlagonal  matrix  corresponding  to  Eq.  (4a), 
set  to  zero  the  epproprlate  elements  of  the  matrix  and 
its  right-hand  aide.  Consider  the  simple  6  *  b  system 


b  c 
a  b  c 


defined  by  Che  minor  grid  indices  k  -  constant  and 
J  ■  1  co  JMAX.  To  achieve  Chls  we  have  used  Che 
following  meChod. 
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If  elemencs  In  rows  3  and  4  are  co  be  blanked,  they  are 
simply  resec  as 
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In  Chis  example,  nondiagonal  elemencs  In  rows  3  and  4 
are  zeroed  so  chac  no  change  is  compuCed  for  (i(/n+1  -  i/*). 
A  similar  creacmenc  must  occur  in  solving  Che  tridiag¬ 
onal  matrix  corresponding  to  Eq.  (4b).  Ultimately,  the 
values  of  v  at  points  3  and  4  are  updated  by  inter¬ 
polation  of  values  from  another  grid. 

Although  the  above  illustration  is  for  a  Laplacian, 
Che  Beam-Warming  (9)  algorithm  for  the  Euler  equations 
is  treated  in  precisely  the  same  manner.  Overset  grid 
results  obtained  using  the  Euler  equation  algorithm 
will  be  given  elsewhere. 

Finally,  the  boundary  condition  must  be  recoded. 
Interpolation  routines  to  update  overset  grid  boundaries 
must  be  provided,  and  modular  coding  is  needed  to 
account  for  the  possibility  that  a  minor  grid  may  have 
different  boundary-condition  treatment  than  the  major 
grid.  For  example,  in  Fig.  2  the  k  =  1  line  in  the 
major  grid  represents  a  periodicity  boundary,  whereas 
the  k  =  1  line  for  the  minor  nose  grid  represents  a 
solid  body  boundary.  Thus,  separate  boundary-condition 
routines  must  usually  be  supplied  for  each  grid,  and 
these  should  be  provisioned  for  in  a  modular  way. 

LOCATING  SPECIAL  GRID  POINTS 

As  part  of  the  overall  grid-generation  package  for 
the  chimera  grid  scheme,  a  program  must  be  included 
that  locates  and  flags  those  points  that  must  be  blanked 
In  the  algorithm  or  that  serve  as  special  boundary 
points.  To  date,  we  have  used  simple  bookkeeping  pro¬ 
cedures  and  allow  the  minor  grids  to  Interact  with  only 
the  major  grid.  The  minor  grids  are  not  allowed  to 
Interact  with  another  body  boundary  (as  shown  in  Fig.  3, 
for  example).  In  time,  these  algorithms  are  expected 
to  be  refined,  generalized,  and  made  more  efficient 
than  current  versions. 

Blanking  Points  Within  a  Body  Boundary 

As  indicated  by  Figs.  1-4,  some  of  the  points  in 
the  main  grid  will  fall  within  the  body  boundary  of  s 
minor  grid.  These  points  must  be  blanked  out.  Because 
the  minor  grid  points  near  the  body  boundary  may  be 
more  finely  spaced  than  the  major  grid  points  that  are 
being  turned  off,  we  should  blank  additional  major  grid 
points  until  the  grid  spaclngs  are  more  compatible. 

Thus,  for  the  situation  shown  in  Fig.  3,  we  may  choose 
to  turn  off  all  major  grid  points  within  some  boundary 


A  boundary  curve  corresponding  co  k  •  constant 
is  defined  within  which  points  are  to  be  blanked.  Call 
chis  curve  C  (see  Fig.  5).  The  origin  to  this  region 
is  defined  by  averaging  all  the  points  on  C.  Ue  then 
search  for  the  point  on  C  which  is  the  farthest  from 
the  origin.  The  distance  from  the  origin  to  this  point 
defines  a  radius,  Rmax.  The  radius  and  angle  for  every 
point  on  curve  C  with  respect  co  the  origin  are  also 
computed  and  stored. 

Points  of  the  major  grid  are  tested  to  determine 
if  they  fall  within  curve  C.  Provision  can  be  made  to 
exclude  points  that  obviously  cannot  be  within  the 
closed-curve  C  by  letting  the  user  input  the  index 
bounds  of  test  points.  This  increases  the  possibility 
for  error,  however,  so  such  an  option  has  not  been  pro¬ 
vided.  The  first  test,  then,  is  to  compute  the  dis¬ 
tance  between  the  major  grid  point  b.-tng  checked  and 
the  origin.  If  this  distance  exceeds  Rmax  then  the 
point  must  lie  outside  of  C;  if  the  point  is  less  than 
Rmax,  its  angle  with  respect  to  the  origin  is  computed 
(see  Fig.  6).  The  points  along  C  are  then  searched 
to  find  that  point  on  C  which  is  closest  in  angle 
measure  to  the  major  grid  point  being  tested.  Points 
on  C  to  either  side  are  then  tested  so  as  to  find 
angle  bounds.  (Because  C  is  closed,  the  angle  will 
Jump  from  -it  to  -Hi  somewhere  between  two  points  of 
C;  if  the  test  point  has  an  angle  within  this  range, 
the  bounding  C-curve  points  are  already  known.)  If 
the  radius  from  the  origin  to  the  test  point  is  less 
than  the  weighted  average  of  the  radii  of  the  bounding 
angles,  the  point  is  within  C.  A  flag  is  then  set  and 
this  information  is  stored.  The  above  test  is  fairly 
reliable  unless  the  body  is  concave  (Fig.  7).  In  this 
case  we  can  partition  C  into  a  C*  and  C"  and  work 
with  two  origins. 

The  previously  described  algorithm  for  locating 
points  within  a  given  boundary  is  only  one  of  many  that 
can  be  envisioned.  An  algorithm  that  can  be  more 
readily  extended  to  three  dimensions  is  sketched, 
although  it  has  not  been  tested  in  computation.  Let 
the  boundary  curve  C  be  represented  by  k  •  constant 
as  before  and  again  assume  that  C  is  continuous  and 
closed.  If  C  corresponds  to  n  *  constant,  and  if 
the  J  index  increases  clockwise  (see  Fig.  8),  then 
represents  an  outward  normal  to  C.  Let  Rp  be 
the  vector  from  the  nearest  point  on  C  to  the  major 
grid  point  being  tested.  If  the  dot  product  of 
with  Rp  is  positive,  the  point  being  tested  is  out¬ 
side  of  the  curve  C;  if  the  dot  product  is  negative, 
the  point  is  inside  of  C. 

location  and  Interpolation  of  Mesh  Overset  Boundary 
Points 


Outer  boundary  data  for  any  minor  grid  must  be 
provided  by  Interpolating  the  solution  of  the  major 
grid.  If  B  is  a  point  on  the  outer  boundary  of  a 
minor  grid,  our  first  task  is  to  locate  the  nearest 
major  grid  point  (or  points)  from  which  Interpolation 
can  take  place.  The  simplest  test  is  to  compute  the 
distance  between  point  B  and  points  of  the  major  grid 
and  to  select  that  major  grid  point  whlrh  is  closest  to 
B.  Again  this  search  can  be  speeded  up  by  excluding 
major  grid  points  a  large  distance  away. 

If  the  boundary  curve  is  smooth,  one  can  also  use 
a  "stencil-walk"  to  speed  the  search.  As  sketched  in 


Fig.  9.  after  the  nearest  point  to  B  is  found,  the 
nearest  point  to  the  next  boundary  point  B'  can  be 
located  by  searching  the  stencil  for  the  point  that  is 
nearest  to  B.  Of  the  stencil  of  nine  points,  find 
that  one  which  Is  nearest  to  8' .  For  this  nearest 
point,  search  its  stencil  for  the  point  closest  to  B* , 
and  so  on,  until  the  point  tested  is  closer  to  B* 
than  any  of  its  stencil  neighbors. 

The  major  grid  will  have  various  "holes”  of 
blanked  points.  The  perimeter  of  these  points  serves 
as  an  inner  boundary  to  the  major  grid,  and  solution 
data  must  be  supplied  to  this  perimeter.  In  this  case 
the  nearest  minor  grid  point  must  be  located  for  inter¬ 
polation  of  minor  grid  data  onto  the  perimeter  of  the 
blanked  major  grid.  A  search  similar  to  that  of  the 
first  approach  outlined  is  used. 

For  Interpolating  the  grid  boundary  data  from  the 
field  of  another  grid  we  have  simply  used  second-order 
Taylor  series  approximations  of  the  form 

*  -  +  A xpx  +  bytyy  +  +  AxAyi^  +  *y y 

(8) 


where 


"  Vt  +  Vn  <etc*> 

and  Ax  and  Ay  are  the  increments  befcven  the  point 
from  which  the  interpolation  is  being  nude  and  the 
point  being  interpolated  for.  A  linear  combination  of 
such  interpolation  formulas  can  be  written  from  nearby 
points  to  ensure  a  very  smooth  interpolation.  We  have 
not  needed  such  a  technique  to  date. 

Because  the  search  of  nearest  points  is  limited 
to  the  index  range  that  excludes  fixed  boundaries,  the 
derivatives  of  and  ip**,  etc.,  can  always  be  formed 
with  central  differencing.  This  simplifies  the  inter¬ 
polation  formula  insofa'r  that  it  does  not  require  the 
U3e  of  special  one-sided  differencing  at  boundaries. 

A  criticism  of  using  such  a  simple  interpolation 
formula  is  that  conservation  is  not  maintained.  That 
la,  although  we  may  difference  the  flow  equations  in 
divergence  form,  the  divergence  property  la  not  numeri¬ 
cally  preserved  if  some  values  are  obtained  by  Inter¬ 
polation.  This  may  be  a  problem  if  a  jump  dlsconClra- 
lty  such  as  a  shock  wave  crosses  the  Interpolation 
boundary.  Perhaps  in  this  case  Interpolated  values  can 
ba  later  adjusted  to  satisfy  a  numerical  line-integral 
of  flux. 

RESULTS 

A a  part  of  a  package  for  generating  overset  grids, 
a  flow-solver  code  has  bean  developed  which  solves  the 
incompressible  stream  function.  Although  such  a  pro¬ 
gram  is  useful  In  its  own  right,  the  application  hare 
is  to  uncover  logic  errors  in  locating  and  flagging 
special  points  in  various  kinds  of  overset  mesh  sys¬ 
tems.  The  ultimate  application  for  this  linear  coda 
will  be  to  use  it  to  assess  the  quality  of  a  given 
generated  overact  mesh  system.  For  example,  the 
Lap lac lan  can  ba  economically  solved  to  help  determine 
whether  an  adequate  masber  of  points  are  provided  at, 
say,  an  airfoil  nose  region  or  whether  an  additional 
grid  should  be  Introduced  to  reaolva  soma  other  feature. 
Specifically,  Eq.  (1)  in  generalized  coordinates  is 


solved  using  tha  algorithm  given  by  Eq.  (4).  Blanked 
points  ara  treated  by  using  the  relaxation  factor  as 
given  by  Eq.  (S) . 

Aa  a  first  teat  of  the  overset  grid  system,  a 
body-fitted  O-grld  waa  superimposed  on  a  stretched 
rectangular  grid  which  served  as  the  major  grid.  Fig¬ 
ure  10a  shows  an  overview  of  an  O-grid  about  an 
NACA  0012  airfoil  overset  on  the  rectangular  grid;  a 
closeup  la  shown  in  Fig.  10b.  The  circular  symbols 
plotted  over  points  of  the  rectangular  grid  indicate 
nearest  major  grid  points  that  are  used  for  interpo¬ 
lating  the  outer  boundary  points  of  the  O-grld.  Here 
the  O-grld  is  considered  to  be  the  minor  grid.  The 
plotted  x  symbols  indicate  points  in  the  rectangular 
grid  that  are  blanked  out.  Although  not  specially 
flagged,  the  perimeter  of  the  blanked-out  points  is 
updated  by  Interpolating  nearby  minor  grid  values.  In 
this  case  all  rectangular  points  within  the  k  »  3 
grid  line  of  the  O-grld  are  blanked.  A  perimeter  of 
rectangular  points  adjacent  to  those  blanked  out  is 
also  flagged.  The  calculated  pressure  distribution  la 
indicated  in  Fig.  10c  for  incompressible  flow  and  la 
shown  compared  with  exact  theory;  the  comparison  la 
satisfactory. 

The  grids  shown  In  Figs.  11a  and  lib  show  a  12Z 
ellipse  modeling  a  nonlifting  biplane  configuration. 

The  lower  boundary  of  the  rectangular  mesh  serves  aa  a 
line  of  symmetry.  Unlike  the  previous  example,  this 
geometry  would  be  poorly  meshed  by  use  of  only  a  single 
O-grld.  This  is  because  the  line  of  symmetry  is  so 
close  that  the  O-grld  would  become  very  distorted.  The 
computed  pressure  coefficient  is  indicated  in  Fig.  11c. 


Another  application  of  the  overset  mesh  system  la 
indicated  by  Figs.  12  and  13  for  an  inlet  with  and 
without  a  centerbody.  In  the  first  case  (Fig.  12),  a 
Stretched  rectangular  mesh  is  used  as  the  major  grid 
and  a  body-fitted  C-grid  is  used  to  resolve  the  inlet. 
The  detailed  grid  view  (Fig.  12a)  shows  blanked-out 
points  and  nearby  interpolation  points  on  the  rectangu¬ 
lar  grid.  Nearby  points  on  the  minor  body-fitted  grid 
which  are  used  to  interpolate  the  perimeter  of  blanked- 
out  rectangular  points  mre  not  shown.  By  choosing 
various  values  of  stream  function  to  be  constant  on  the 
body  boundary  we  can  control  the  mass  flow  into  tha 
inlet.  Computed  streamlines  are  shown  In  Figs.  12b 
and  12c  for  mass  flow  rates  that  cause  spillage 
(Fig.  12b)  and  suctloi  Fig.  12c).  An  additional  aat 
of  grida  and  computational  results  are  shown  In 
Figs.  13a-13c.  In  this  latter  case  the  inlet  has  a 
centerbody  that  is  fitted  by  shearing  Che  major  rec¬ 
tangular  grid  over  this  boundary. 


A  final  example  shows  an  airfoil  and  flap  arrange¬ 
ment  using  0-grlds  for  both  the  major  and  minor  grids. 
This  arrangement.  In  which  the  flap  It  tucked  In  below 
tha  main  airfoil,  esnnot  be  nicely  fitted  by  s  single 
mapping  with  cuts.  Overviews  of  the  grids  end  computed 
results  ar#. Indicated  In  Figs.  14a-14d.  The  circula¬ 
tion  of  each  profile  is  found  automatically  In  these 
cases  by  requiring  a  match  of  tralllng-edga  praasuras 
aa  a  Kucta  condition.  Using  the  Bernoulli  equation, 
tha  condition  of  setting  the  preesures  at  the  point 
above  the  trailing  edge,  u,  to  the  point  below,  1, 
gives  tha  raise  ion 


Ko*  ♦ 


•t1)**  1 

y  n  upper 


f(nj  ♦ 


n1 lower 


(9) 


This  relation  uaaa  tha  lnvlscld  txngency-bouodary  con¬ 
dition  that  la  constant  on  the  body  ao  “  0;  1  and  n 


x 
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represent  the  curvilinear  body-conforming  coordinates 
with  (  circumferential  and  n  directed  outward. 
Differencing  this  relation  at  the  trailing  edge  results 
in  a  solution  for  the  stream  function  on  the  body  in 
terms  of  neighboring  ^-values  at  the  trailing  edge. 

CONCLUDING  REMARKS 

It  has  been  known  for  a  long  time  that  multiple- 
grid  systems  work  and  that  they  offer  computational 
advantages.  Although  using  multiple  grid  systems  Is 
conceptually  straightforward,  many  experienced  prac¬ 
titioners  of  computational  fluid  dynamics  (CFD)  who 
have  coded  such  programs  would  likely  agree  that  all 
of  the  extra  bookkeeping,  although  simple,  tends  to  be 
tedious.  It  Is  clear,  however,  that  some  type  of  mul¬ 
tiple  grid  will  be  utilized.  A  major  task  in  CFD, 
then,  will  be  to  design  computer  codes  that  can  use 
multiple  grids  and  still  remain  usable  and  readable. 

In  this  paper  we  have  described  a  study  of  a 
chimera  grid  code  architecture  using  overset  grids. 

We  are  finding  chat  such  a  grid  system  need  not  overly 
complicate  the  computer  code,  that  it  can  be  handled 
within  a  general  framework,  and  that  It  can  simplify 
the  Cask  of  grid  generation. 
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Fig.  1  Overset  grids  for  airfoil  with  flap  Fig.  2  Overset  grids  for  cascade  blades  with  minor 

grid  used  Co  resolve  nose  of  blade 


Fig.  3  C-grld  for  cascade,  showing  extensive  skswneaa 


Fig.  4  Airfoil  with  leading-  and  trailing-edge  flaps 
resolved  with  overset  grids;  note  that  the  grid  spac- 
ings  are  not  shown  to  correct  scale 


rig.  5  Curve  C  which  defines  region  of  blanked-out 
major  grid  points 


Fig.  6  Test  procedure  for  locating  major  grid  points 
that  lie  within  curve  C 


Fig.  7  Partition  of  curve  C  for  case  of  concave 
boundary 
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a)  Overview 


c 


Pressure  distribution  compared  with  exact  linear  solution 


b) 


Detailed  view 


Fig.  10  0-grld  about  NACA  0012  airfoil  on  Cartesian  major  grid. 
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I . _ Abstrjv:  t_ 

An  implicit,  approximate-factorization,  finite- 
difference  algorithm  has  beer,  developed  for  the  com¬ 
putation  of  unsteady,  inviacid  transonic  flows  in 
two  and  three  dimensions.  The  computer  program 
solves  the  full-potencial  equation  in  generalized 
coordinates  In  conservation- law  form  in  order  to 
properly  capture  shock -wave  position  and  speed.  A 
body-fitted  coordinate  system  is  employed  for  the 
simple  and  accurate  treatment  of  boundary  conditions 
on  the  body  surface.  The  time-accurate  algorithm  is 
modified  to  a  conventional  ADI  relaxation  scheme  f- 
steady-state  computations.  Results  from  two-  and 
three-dimensional  steady  and  two-dimensional  unsteady 
calculations  are  compared  with  existing  methods . 


II . _ Introduction 

Modern  transport,  military  and  rotor  aircraft 
routinely  operate  in  the  transonic  flight  regime. 

In  transonic  flight,  the  aerodynamic  forces  are 
extremely  sensitive  to  small  perturbations  in  the 
motion  of  the  vehicle.  Unsteady  flow  adjustment  is 
very  slow  as  upstream  propagation  of  Information  is 
restricted  by  locally  high  subsonic  or  lupersonic 
flow  resulting  in  large  phase  lags.  In  addition, 
the  existence  of  embedded  6hock  waves  with  the  poten¬ 
tial  for  large  excursions  further  complicates  the 
transonic-flow  environment.  Thus,  it  is  not  sur¬ 
prising  that  the  most  critical  aeroelastic  phenomena 
occur  In  this  flight  regime.  And  while  the  super¬ 
critical  airfoil  sections  are  more  aerodynamical ly 
efficient  than  c ■  iventicnal.  airfoils,  they  are  more 
susceptible  ,c  flutter.1  Linearized,  unsteady  sub 
sonic  aerodynamic  theory  Is  incapable  of  predicting 
these  complicated  flows.  In  addition,  transonic 
prediction  methods  that  use  the  harmonic  approach  to 
compute  a  small  perturbation  from  nonlinear,  mean 
steady-state  flow  have  limited  applicability.  They 
are  valid  only  for  smal 1- amplitude,  high-frequency 
flows  and  cannot:  for  shock-wave  motions. 

Thus,  it  is  of  paramount,  importance  to  develop 
methods  that  can  properly  compute  shock  waves  and 
their  motions  if  their  /ole  in  transonic,  aeroelas 
tic  stability  is  to  be  accounted  for.2 

Unsteady,  transonic,  aerodynamic  method,.,  ba-  d 
on  various  run linear  I ! -d i nt  vbance  potential 
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equations  have  been  developed  by  several  research¬ 
ers.  Ballhaus  and  Steger3  and  Ballhaus  and 
Goorjian*'  have  developed  an  efficient  method  for 
solving  the  two-dimensional  low-frequency  transonic 
small-disturbance  equation.  The  resultant  code, 
LTRAN2,  has  been  extensively  used  in  spite  of  its 
various  limitations.  Recently,  methods  that  moder¬ 
ately  extend  the  frequency  and  Mach  number  range  of 
this  equation  have  been  developed.5-7  In  addition, 
the  effect  of  viscous  corrections  has  been  examined 
by  Rizzetta.8  In  another  effort,  Borland,  Ri2zetta, 
and  Yoshihara9  have  recently  reported  on  an  algo¬ 
rithm  for  a  modified  form  of  the  small-disturbance 
equation  in  three  dimensions  with  no  frequency  limi¬ 
tations.  This  code  is  currently  being  used  to  com¬ 
pute  unsteady  aerodynamics  for  transonic  flutter 
analyses.10*  1 

A  practical  tool  for  the  computation  of 
unsteady,  invlscid  transonic  flow  about  complex 
configurations  must  accurately  simulate  the  signifi¬ 
cant  physics  and  be  computat ionally  efficient  for 
routine  use  in  engineering  analysis.  Unsteady, 
full-potential  theory  can  satisfactorily  replace 
Euler  equation  solutions  if  the  shock  waves  are 
sufficiently  weak  and  yet  maintain  computer  time 
and  storage  requirements  similar  to  the  simplified 
small-disturbance  theory.  The  unsteady,  full- 
potential  equation  has  been  solved  using  fully 
Implicit  methods  by  Goorjlan,12  Steger  and 
Caradonna.13  Chipman  and  Jameson, *  *  and  Sankar 
et  al.15*  8  The  Chipman  and  Jameson  procedure 
solves  a  system  of  two  equations  for  the  two 
unknowns,  density  and  velocity  potential,  using  an 
approximate-factorization  scheme.  The  method  of 
Goorjian  uses  an  ADI  scheme  to  solve  a  scalar  equa¬ 
tion  for  the  velocity  potential  by  employing  a 
time-linearization  of  the  density  similar  to  the 
present  approach.  To  date  both  of  these  methods 
have  only  been  applied  to  two-dimensional  nonlifting 
pulsating  airfoils  using  a  simple,  sheared  mesh. 

The  method  of  Sankar  et  al.  uses  the  strongly 
implicit  procedure  and  has  been  applied  to  compute 
steady  and  unsteady  flow  over  wings  using  a  sheared, 
parabolic  coordinate  sysiem. 

The  present  procedure  is  an  extension  of  the 
Steger  and  Caradonna  work13  which  used  a  stretched 
r3rte°.ian  grid  with  small-disturbance  planar  bound¬ 
ary  conditions.  Currently,  a  generalized  body- 
fitted  coordinate  system  is  employed  that  can  be 
adapted  to  wings,  wing-body  combinations,  and 
bodies  of  revolution.  This  method  is  unique  as  it 
can  use  the  mesh-point-efficient  spherical  grids 
that  are  presently  being  developed  in  conjunction 
with  this  code.  Steger  and  Caradonna  reported  that 
the  unsteady  scheme  was  a  very  poor  relaxation 
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algorithm.  In  thla  work,  tha  achaaa  la  nod If lad  to 
an  ADI  relaxation  procedure  for  computing  faat, 
steady-state  aolutlona  by  a  al^la  flag. 

The  governing  a qua t Iona  and  boundary  condi- 
tlona  are  dlacuaaad  In  Secdon  II.  In  Section  III 
tha  conaervatlve  differencing  and  numerical  algo¬ 
rithm  for  the  aolutlon  of  tha  flnlta-dlfferanca 
equation#  are  praaantad.  In  Section  IV,  raaulta  are 
prea anted  for  varlou*  two-  and  three-dlaenalonal 
caaaa  to  verify  the  algorithm.  Finally,  concluding 
remarks  are  Bade  In  Section  V. 


III.  Mathematical  Formulation 


Full-Potential  Equation 


with  tha  Bernoulli  relation,  Eq.  (2),  trana forming 
aa 

P  -  |l  +  [»!  -  2»t  -  (U  +  ct)F? 

-  (V  ♦  ne)#n  -  (W  ♦  <t>Ajj1/T_1 


0  -  et  +  A^t  +  A^n  +  A,AC 
v  -  nt  -*•  *  K,*rt  +  At*t 
W  -  tt  t  a5a?  +  t  Asac 


The  three-dlaenalonal,  unateady,  full-potential 
equation  In  atrong  conservation-law  form  la  given  by 

It +  h  (<**> +  <pV  +  h  (p4*)  ’ 0  (1> 

where  the  density  p  la  determined  from  the 
unsteady  Bernoulli  relation,  for  steady,  uniform 
Incoming  flow 

P  -  fl  +  *-=-5-  -  2*t  -  A2  -  ♦*  -  A2)}1''*'1  (2) 


Ai  "  +  «y  +  <z 

A*  *  nx  +  ny  +  nz 

A>  "  <  +  ^  +  <z 

A„  -  ?2n2  +  C2n2  +  52n2 

*  x  x  y  y  i  * 

A.  “  eV  +  c2t2  +  t2;2 
5  <  y  y  *  z 

A,  -  n2;2  +  n2C2  +  n2;2 

•  x*x  y  y  a** 


The  density  p  and  velocity  components  Ax>  Ay, 
and  A2  are  nondlmenslonsllsed  by  the  free-sereaa 
p„  and  speed  of  sound  a„.  The  Cartesian  coordi¬ 
nates  x,  y,  and  z  are  nondlmanslonalized  by  a 
reference  length  t,  such  as  the  airfoil  chord  c, 
and  the  time  t  la  referenced  to  a„/l. 

Equations  (1)  and  (2)  express  mass  conservation 
for  unsteady,  lnvlacld,  lsentroplc,  and  lrrotatlonal 
flow.  The  corresponding  shock- jump  conditions  can 
be  a  suitable  approximation  to  the  Kanklne-Hugonlot 
relations  for  many  transonic  flow  applications  If 
the  shock  waves  are  sufficiently  weak.  In  deriving 
Eq.  (2)  Che  far-fleld  la  assumed  to  be  steady. 

Transformation  of  the  Equations 

The  treatment  of  arbitrary  body  boundaries  Is 
usually  made  more  convenient  by  the  use  of  a  coor¬ 
dinate  transformation  which  maps  the  body  surface  to 
a  rectangular  coordinate  surface  in  the  transformed 
plana.  Boundary  conditions  at  the  body  surface  can 
then  be  simply  and  accurately  treated.  In  addition, 
these  mappings  can  be  used  to  cluster  grid  points  In 
regions  of  the  flow  with  high  gradients,  thereby 
enhancing  mserlcal  solution  accuracy.  A  general 
Independent  variable  transformation  is  Indicated  by 

K  -  5(x,y.z.t) 

n  “  nCx.y.z.t) 


J  "  “  n.C„)  +  n  ({c  -  Sc) 

*  y  *  *  y  xzy  y  z 

♦  StV* '  W  (6lr) 

Here,  U,  V,  and  W  are  contravarlant  velocities 
along  the  t,  n.  and  ;  directions,  respectively, 

Aj  -  A,  are  metric  quantities,  and  J  Is  the 
Jacobian  of  the  transformation  I 3(t,n,c)/3(x,y,r) I 
The  metric  quantities  In  Eqs.  (6a)  and  (6b)  are 
evaluated  using  the  following  metric  Identities 


’  J(xnV  xtV 

nx  "  J(5,c*c  ■  Vt’ 

"y  "  J(xsWs) 
nz  “  J <xcyt  ’  W 


;X  *  -,(yc*n-ynV 

■y  '  J(Vi  '  Vn> 

(.  •  J  (xry  -  x  y.  ) 

L  »  -x  r  -  y  T  -  zC 
St  T  X  7  T  Sy  T  z 

*  -  x  n  -  y  n  -  t  rt 
t  1  X  T  y  T  z 


Boundary  Conditions 

At  the  body  surface,  flow  tangency  Is  required 
il.e.,  no  flow  through  the  body)  This  Imposed  by 
setting  the  contravarlant  veloeltv  In  the  c,  direc 
tlon  to  zero,  viz 


C  -  (fx,y,z,t) 


W  -  t  +  Ast,  +  A.An  4  A,*.  -  0 


The  strong  conservation- law  form  of  Eq.  (1)  is  main¬ 
tained  by  axprasslng  it  as 


A  similar  condition  (i.e.,  V  -  <■  f  t*  Impoaed  at  the 
symmetry  plane  tor  win*  cab'h. 

For  lifting  canes  ,  a  1  low  ant.  «•  must  be  made  roi 

a  Jump  of  potential  across  a  wake  like  cot.  In 
unsteady  flow,  vortlcltv  le  continuously  abed  from 
the  wing.  In  a  potential  formulate  cm  this  ib 
approximately  modeled  with  a  cut  aligned  with  the 


(12) 


shear  layer.  Here  we  further  Idealize  this  flow  by 
keeping  the  cut  In  the  nean  chord  line  plane.  By 
laposlng  the  usual  shear-layer  assumptions >  an 
equation  for  the  Juap  In  potential. 


r  -  -  *, 


along  the  cut  can  be  derived.  Across  an  lnvlscld 
shear  layer,  normal  velocity  and  pressure  are  con¬ 
tinuous.  Since  lsentroplc  flow  has  been  assumed, 
the  density  It  also  taken  to  be  continuous.  The 
Bernoulli  relation  together  with  the  continuity  of 
density  requites  that 

r  +  <v>r  +  <w>r  -  o  (9 

T  ri  4 

where  (V)  and  (W)  are  the  averages  of  the  contra- 
variant  velocities  above  and  below  the  wake,  i.e.. 


<V>  -  j  <Vu  +  Vt> 


<„>  "  J  <W„  +  V 


This  convection  equation  for  r  is  Imposed  in  the 
wake. 

At  distances  far  from  the  body  the  flow  is 
required  to  be  free  stream,  which,  with  the  use  of 
p„  and  a„  as  the  references,  is 

♦  “ 


+  V3n  +  wa;j 

Finally,  to  better  facilitate  the  application 
of  approximate  factorization  techniques,  the  cross- 
derivative  terms  are  lagged  In  time  by  rewriting  the 
Implicit  spatial  derivatives  In  the  form 

^(pU)”*1  -  35(pAl)n35<*n+1  -  ♦“)  +  35(pU)n  +  0(h)| 
3  (pV)“+l  -  3  (pAa)n3  (dn+1  -  ♦“)  +  3  ( pV)“  +  0(h)! 


3^ (pW)n"*1  -  3c(pAJ)n3|.($n+l  -  *n)  +  3c(pW)“  +  0(h)J 


Combining  Eqs.  (10)-(13)  yields  the  conservative 
time-dlscretlzed  form  of  the  full-potential  equation 
involving  only  the  potential  at  the  new  time  level 

jl  +  h(0n3?  +  V1^  +  Wn3?)  -  |i  l3?(pAl)n3t 

+  3n(£A2)n5ii  +  3;(^A3)I>3|.j|(tn+1  -  *n) 


(*n  -  ♦"  *)  +  2-—  (*n  -  2<pn-1  +  *n-2) 
6 


+  —  (p“  -  Pn_1)  +  h  ^ -  (l^'V  +  Vn~l9 

Bn  Bn  E  n 


where  M„  is  the  free-stream  Mach  number.  For 
steady  lifting  cases,  the  velocity  potential  at  the 
outer  boundary  is  updated  with  the  usual  compressi¬ 
ble  vortex  solution  with  strength  T. 


Ill.  Numerical  Algorithm 


Temporal  Differencing 

A  first-order-time-accurate  approximation  to 
Eq.  (3)  is  obtained  using  Euler  backward  differencing 

S"+I  -  p"  +  h(3.(p(l)n+l  +  3  (3V)n+1  +  3  (pW)n+1]  -  0 


where  h  -  At  and  p  *  p/J.  To  solve  a  scalar  sys¬ 
tem  of  equations  for  the  potential  at  each  new  time 
level  the  density  is  linearized  about  old  time 
levels.  The  density  coefficients  in  the  spatial 
derivative  terms  are  lagged  to  level  n.  A  linear¬ 
ization  resulting  in  a  conservative  differencing  of 
the  time  derivative  ot  the  density  is  obtained  by 
noting  that  p  »  p($)  and  using  a  Taylor  series 
expansion 

-n+i  -n  f  n  (\  3p\n.  n+i  ,n,T 

p  -  (i  *  [P  t  ^  (*  '  ♦  >J 


where  3p/3li  is  a  differential  operator  obtained 
from  the  Bernoulli  relation 


+  Wn'13c)($"  -  4>n_1 )  +  [3t(pU)n  +  3n(pV)n 

+  3;(pW)n]  (14) 

where  i  -  p2wy/J. 

Spatial  Differencing 

An  optionally  first-  or  second-order-accurate 
spatial  differencing  of  Eq.  (14)  is  given  by  (note: 
the  spatial  indices  are  suppressed  here  for 
convenience) 

li  +  h(uV  +  vn«  +  wn«  )  -  ^  (pAl)n3 
(  t  n  C  ;n  C  1  C 


nv)<*r 


+  6r)(pAJ)"6ri  +  fi^tpAj) 


(♦“  -  t"'1)  +  (♦"  -  2*"'1  +  *n-2) 


+  -  <p"  -  Bn")  +  h  SJ-  |un-’4f  ♦  Vn-*6 

b"  8"  5  n 


♦  v"-\l(*n  -  ♦"-*)  +  [?t(f )"  +  4n(pV)n 

+  3  ^  ( pC ) n J  Ob 


where 


-  Vi/ A/,(^)<4  -  ^>J 

*n<PA2>!n*k  -  )(*k+i  "  V 

-  -  Vx>] 

3c(pAJ)?c4t  -  |AJjt+l/2Jj|l/2(  4  ')(*t+i  ■  ♦l) 

-  A>i-l/2JI-l/^f^^H<♦l  -  wl 


(16) 


and 


“e(V)  "  Jj+i/a4+1/J°J+i/2 

JJ-i/j4-i/2^J-i/2 

Vp*>  ■  J^/fkt‘rk)v./i 

'  ^-l/2(Ck  Vk73)Vi/a 

V£S)  • 


_  j-i  /p*  *  °t-Ae 

l-l/2\  2  }  1-1/3 


(17) 


The  contravariant  velocities  at  the  meBh  half-points 
are  determined  from 


“j  +  l/2  ’  ?tj  +  l/2  +  Alj  +  l/2(*j  +  l  “  V 


+ 1  [A*j+.Vj+i +  N4-.#j] 
+  1  [A*J+lVj+i  +  A*jVj 


with  similar  treatment  for  Dj-i/2.  Vk-i/i .  and 
fit-! /j.  The  metric  terma  at  these  points  are 
obtained  using  simple  everagea  (e.g. , 

AM+t/»  ■  (Axj  +  Aij+l)/2).  Here  A?  ■  An  "  AC  “  1 
and  only  the  varying  Indices  are  indicated.  The 
derivative  terms  in  the  contravariant  velocities 
expressions  are  replaced  with  the  second -order- 
accurate  central  difference  approximations 

1 


♦r|  *M 

'j 


ci4  ~  -rj  7  (4+i "  4j-i 5 


♦n  4  Vk  "  I  (*k+»  _  ♦k-i) 


♦c  .  4  V»  ■  7  (*1+1 


♦t-i} 


(19) 


Stability  is  maintained  in  supersonic  regions 
by  introducing  an  artificial  viscosity  term  through 
the  use  of  an  upwind  bias  of  the  density 


PJ+l/2 


(1  -  /->(P^2~"j) 


+  v 


J+ 


i+ih‘ 

i/2[  2 


pJ+r  +  "S-51  pJ-i+>r] 

(20) 


where  r  *  0  or  1  for  U  <  or  >  0.  The  parameter 
8  •  2  for  second-order  spatial  accuracy  in  super¬ 
sonic  regions,  and  0-0  for  first-order  accuracy. 
The  switching  paraaMter  v  is  defined  by 


j max[ 1  -  (p/p*)j,0)C 


U.^  ,  >0 

j+1/2 


j+l/l 


j  max [ 1  -  (p/p*)j+1,0]C  U 


V>/2 


<  0 


with  1  <  C  <  10.  Note  that  upwindlng  is  used  in 
the  t  direction  only  at  present.  For  the  compu¬ 
tations  performed  to  date  this  has  been  satisfac¬ 
tory,  but  in  general  it  will  be  necessary  to  add  it 
in  the  other  directions  as  well  (see  kef.  17  for 
the  extension). 

The  metric  quantities  on  the  right-hand  side  of 
Eq.  (17)  are  computed  using  three-point,  second- 
order-accurate  difference  expressions 


"C  7  lV»  •  xJ-.) 


\  -  2  (Vm  -  V.’ 


(21) 


Vi/2 

“  ntk+i/ 

2  +  A*k+i/r(*k+i  "  V  I 

+  }  [a 

»k+i4E*k+i  +  A|*k4C4k] 

+  7  [A 

‘k+i*i*k+i  +  A*kMk] 

fil+i/2 

“  Ct*+i/ 

2  +  Asi+i/2^l+i  “  *0 

+  7  [A 

’i+i4E*i+i  +  Asl44*t] 

+  7  [* 

*t+l4n*t+l  +  A‘t4n*tJ 

xc  *  7  <\+.  -  xt-.)  ) 

(18)  with  similar  treatment  for  the  y  and  i  metric 

terms.  At  the  boundaries,  3-point,  one-sided  dif¬ 
ference  expressions  are  used. 

It  is  necessary  to  subtract  a  numerical  trunca¬ 
tion  error  term  because  of  an  incomplete  metric 
cancel  1st  ion . 1 5 ’ 1 *• 1 ’  Tf  the  velocity  and  density 
are  set  to  free-stream  conditions  then  the  terms  in 
Eq.  (17)  will  not  be  identically  zero,  but  will  be 
proportional  to  the  numerical  truncation  error  in 
the  differencing  of  the  metrics.  This  is  caused  by 
the  choice  of  differencing  used  for  the  metric  quan¬ 
tities  in  Eq.  (21)  and  the  spatial  derivative  terms 
in  Eq.  (17).  The  error  term  can  be  quite  appreciable 
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for  a  highly  atretchad  grid.  Let  R«,  repraaent  the 
truncation  error  term  obtained  by  aettlng  4  “  H_x 
and  p  -  1  initially  and  computing  the  terma  in 
Eq.  (17).  i.e., 

*-  ■  »«(¥) ♦  k(¥)  •  \(¥) 

This  error  term  la  then  subtracted  from  the  right- 
hand  side  of  Eq.  (IS)  at  each  time  step.  This 
correction  is  particularly  effective  in  the  highly 
stretched,  coarse-grid  far-fleld  where  the  solution 
is  close  to  free  stream.  Near  the  airfoil  the  grid 
resolution  should  be  sufficiently  fine  that  R„  Is 
negligible. 

Approximate  Factorization 

To  avoid  costly  matrix  Inversions  at  each  time 
level  Eq.  (15)  is  approximately  factored  into 
L^,  Ln,  and  Lj.  operators 


x  [i  +  hV*V  -  6  (pA2) 

L  n  gn  n 


|n  «5(PAi>n^j 

—  6  (pA,)nS  ] 
gn  n  2  „J 

^  6r(pA,)ni  ( 

i"  5  3 


1  +  hW\  '  Tn  «5<PA3)n2;J(4n+l 


(♦"  -  4n_1)  +  (4n  -  2*n'1  +  4n_2) 

6° 


+  —  (p"  -  p"~l)  +  h  (Un_  1 6  +  Vn-l6 

gn  gn  E  n 

P 

+  ?c(pW)n  -  R^J 

This  equation  has  the  form 


apace-tlme  derivative  terma  on  the  left-hand  side 
with  a  simple  flag  giving 

{I  -  hJe(pA1)“Jel[I  -  h3n(pAj)n?n](l  -  h5c(pA,)nSt) 

-  (4n+1  -  ♦“>  -  *>[*?(^)“  +  ypv>”  +  ?c(pe,)n  -  *„] 


This  is  satisfactory  for  subcritical  and  slightly 
critical  cases,  but  additional  temporal  damping  will 
be  required  for  cases  involving  large  regions  of 
supersonic  flow.  For  these  cases,  the  space-time 
derivatives  on  the  left-hand  side  of  Eq.  (23)  may  be 
turned  back  on  with  proper  upwind  bias  to  provide 
the  necessary  temporal  damping  and  yield  steady 
state  performance  that  approaches  the  AF2  scheme. 

Boundary  Condition  Implementation 

The  body-surface  flow  tangency  condition  is 
imposed  by  using  Eq.  (8)  to  derive  a  space-time 
extrapolation  procedure  for  updating  the  solution  on 
the  body.  A  second-order  accurate,  finite- 
difference  approximation  to  Eq.  (8)  can  be  derived 
by  replacing  the  (  and  q  spatial  derivative  terms 
with  central  difference  operators  and  the  t  deriv¬ 
ative  with  a  three-point  one-sided  difference 
expression.  The  resulting  equation  is  written  ln 
delta  form  and  approximately  factored  to  yield 


(*  -  if;  ‘t)  (■  -  ifeh  »: 


I 

3  '^2 


♦r> + 


3X7  [;t  +  <MC'  +  MnH2+1l 


where  2.  At  the  trailing  edge, 

two-point  one-sided  difference  operators  are  used 
for  6j.  Similarly,  at  the  symmetry  plane  a  one¬ 
sided  operator  is  used  for  6n.  After  the  interior 
flow  field  has  been  updated,  Eq.  (27)  is  then  solved 
to  update  the  solution  on  the  body. 


L?LnL;(4  -  ♦")  -  R 

and  it  is  implemented  as  algorithm  as 
L^A4*  »  R  | 


L 


Lj.44n+1  -  *♦**  I 

,n+i  ,n  ,  .  ,n  I 

•4  »  4  +  / 

The  algorithm  Eq.  (2S)  requires  only  a  series 
of  scalar,  trldlagonal  inversions  and  it  is  there¬ 
fore  very  efficiently  solved.  Computer  storage  for 
three  levels  of  4  and  one  level  of  p  are 
required. 

Steady  State  Algorithm 

A  fast,  steady  state  ADI  relaxation  algorithm 
is  readily  obtained  from  Eq.  (23)  by  turning  off  the 
unateady  tanas  on  the  right-hand  side  and  the 


The  governing  equation  is  solved  on  the  sym¬ 
metry  plane  for  wing  cases  using  the  warped,  cylin¬ 
drical  coordinate  system  shown  in  Fig.  1.  Along 
this  boundary  flow  tangency  is  imposed  by  setting 
the  fluxes  on  each  side  of  the  wall  (k  ■  1)  to  be 
equal  and  opposite,  viz 


(Si)  .  _(£V\ 

'k-l/2  'J'k-3/2 


For  this  grid  topology,  the  solution  for  the  veloc¬ 
ity  on  the  wing  extension  (i.e.,  the  flat-plate  sec¬ 
tion  beyond  the  wing  tip)  is  determined  by  averaging 
the  solution  from  known  updated  values  one  surface 
away  (£  •  2).  In  particular,  the  potential  is  set 
to  be  the  free-stream  value  on  the  flat-plate  sec¬ 
tion  (i.e.,  4„(x,0,z))  plus  the  average  of  the  per¬ 
turbations  from  free-stream  values  at  the  points 
Immediately  above  and  below.  The  density  is  also 
extrapolated  and  averaged  on  the  wing  extension. 

The  new  values  of  the  circulation  r  ln  the 
wake  are  obtained  by  solving  Eq.  (9)  after  the 
velocity  potential  has  been  updated  in  the  flow 
field  and  on  the  body  surface.  Equation  (9)  ia 


solved  using  a  finite-difference  approximation  where 
3n  and  »|j  ara  replaced  by  tbs  central  and  backward 
difference  operators!  t.  and  respectively.  Tbs 
aquation  la  written  In  delta  fora  as 

(I  +  At<w)a  ♦  At (V>°«  )Ar“ 
n 


-  -At«V)n«  +  <W>nf„ )r“  +  At<W)aAr?  ,  (28) 

n  c  i-i 

Equation  (28)  Is  used  to  update  T  to  the  new  tlaa 
level.  The  tralllng-edge  value  of  T  Is  obtained 
by  setting  it  equal  to  the  new  value  of  the  jmp  in 
potential  there.  The  solution  to  Eq.  (28)  Is  than 
coaputed  by  marching  in  the  i  direction,  thus 
solving  tridlagonals  in  the  n  direction  at  each 
C  -  constant  line. 

For  nonlifting  flows,  the  outer  boundary 
remains  fixed  at  free-streaa  conditions.  For  steady 
lifting  flows,  the  values  of  the  potential  are 
updated  using  the  compressible  vortex  solution 

♦  -  M-*  +  h  6  (29) 

In  the  preliminary  calculations  obtained,  the  outer 
boundary  values  for  the  potential  in  the  unsteady 
case  remain  unchanged  from  their  steady-state 
values.  Care  is  taken  to  place  the  boundary  suffi¬ 
ciently  far  away  to  prevent  reflected  waves  from 
contaminating  the  solution  at  the  body. 


Density  Update 

After  the  velocity  potential  has  been  updated 
in  the  entire  flow  field  using  Eq.  (25)  and  on  the 
boundaries  and  wake  using  Eqs.  (27),  (28),  and  (29), 
the  density  is  updated.  The  density  is  computed 
from  the  Bernoulli  relation,  Eq.  (9),  where  the 
spatial  derivatives .are  replaced  with  the  difference 
expressions  of  Eq.  (19)  and  the  time  derivative  with 
a  two-point  backward  difference  operator.  At  the 
body  surface,  an  expression  for  the  (  derivative 
is  determined  by  solving  Eq.  (9)  for  i.e. , 
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and  is  also  used  in  the  computation  of  the  contra- 
variant  velocities  on  the  body.  At  the  wall  boundary 
(y  •  0),  a  two-point,  one-sided  difference  expression 
was  used  for  the  n  derivative.  All  exponential 
functions  were  eliminated  by  using  binomial 
expansions. 


Grid  Generation 


A  nusber  of  grid-generation  programs  have  been 
used  in  the  present  investigation  because  of  the 
wide  range  of  geometries  considered.  The  0-type 
grids  used  for  the  airfoil  calculations  were 
obtained  using  the  grid-generation  program  GRAPE 
(grids  about  airfoils  using  Poisson's  equation).20 
This  grid-generation  schme  maerlcally  generates 
solutions  to  Poisson's  equation  to  establish  regular 
and  smooth  finite-difference  meshes  around  arbitrary 
two-dimensional  bodies.  The  inhomogeneous  terms  are 
automatically  chosen  to  control  the  mesh  point  spac¬ 
ing  adjacent  to  the  boundaries  and  the  angles  with 
which  mesh  lines  Intersect  the  boundaries.  The 
aquations  are  transformed  to  and  solved  in  the 


computational  domain  using  successive  overrelaxation 
(SIDE) . 

The  grid  genaratlon  program  used  for  wing  cases 
that  use  the  warped  cylindrical  coordinates  was 
presented  in  Ref .  17.  The  finite-difference  mesh 
is  generated  by  employing  a  standard  two-dimensional 
grid-generation  schasM  similar  to  GEAPE.  Mtaerlcal 
solutions  to  the  elliptic,  partial  differential 
grid-generation  equations  are  obtained  in  each 
a panel a a  plane  used  as  a  defining  station  for  the 
wing.  The  equations  are  transformed  to  and  solved 
in  the  computational  domain  using  a  fast 
approximate-factorization  algorithm.  This  estab¬ 
lishes  values  for  x  and  z  in  each  spanwise  plane. 
The  coordinate  values  in  the  spanwise  direction 
(y  values)  are  computed  from  a  stretching  formula 
that  in  its  simplest  form  gives  equal  spacing  over 
the  wing  with  relatively  rapid  stretching  beyond 
the  tip.  For  the  wing  extension,  a  flat-plate 
section  is  used. 

For  the  three-dimensional  flows  using  warped 
spherical  coordinates,  a  hyperbolic  grid-generation 
procedure  was  used  (Steger,  J.  L. ,  Jesperaen,  D. , 
and  Strlgberger,  J.,  private  communication).  This 
procedure  is  a  three-dimensional  extension  of  the 
scheme  devised  by  Steger  and  Chaussee21  in  which  a 
system  of  hyperbolic  equations  is  solved  using  an 
implicit,  marching  finite-difference  scheme.  Two 
of  the  equations  are  derived  from  orthogonality 
conditions  and  the  third  relation  i<-  a  specifica¬ 
tion  of  the  mesh  cell  volume.  Since  this  Is  a 
noniterative  algorithm,  very  fast  grid  generation 
is  obtained. 


IV.  Results 


The  algorithm  Eq.  25)  has  been  coded  Into  a 
computer  program  named  TUNA  (transonic  unsteady 
aerodynamics).  The  three-dimensional  computer  code 
functions  as  a  two-dimensional  code  by  simply 
setting  8  flag.  Consequently,  steady-state  calcu¬ 
lations  in  two  dimensions  were  first  performed  to 
verify  the  accuracy  of  the  steady  and  unsteady 
algorithms.  The  solutions  obtained  with  the  present 
method,  for  several  standard  test  cases,  have  been 
compared  with  the  2D  steady,  transonic  computer 
code  TAIR  (transonic  airfoil  analysis) 2 2 ‘ 22  which 
solves  the  steady  full-potential  equation  using  the 
approximate  factorization  scheme  AF2. 

A  comparison  of  solutions  for  a  subcrltlcal 
nonlifting  test  case  Is  shown  in  Fig.  2.  The  pres¬ 
sure  distribution  on  the  upper  surface  of  an  NACA 
0012  airfoil  at  *  0.72  and  o  ■  0°  shows 
excellent  agreement  between  the  two  codes.  Both 
codes  used  100«31  0-mesh  type  grids  with  TAIR  being 
Internally  generated.  The  grid  used  In  the  present 
method  was  generated  using  GRAPE  (Fig.  3).  This 
steady-state  computation  was  performed  using  the 
steady  and  unsteady  algorithms.  An  optimum  set  of 
acceleration  parameters  was  found  for  the  steady 
algorithm  and  an  optimum  time  step  for  the  unsteady 
algorithm.  The  residual  histories  for  these  compu¬ 
tations  are  shown  in  Fig.  4.  The  new  steady  scheme 
displays  rapid  steady-state  convergence  and  la  at 
least  an  order  of  magnitude  faster  than  the  unsteady 
scheme.  The  present  method  dropped  the  maximum 
residual  four  orders  of  magnitude  in  43  iterations 
compared  to  74  for  TAIR.  Past  experience  has  shown 
that  the  ADI  convergence  rate  Is  about  twice  as 
fast  as  AF2  for  subcrltlcal  flows. 
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The  aubcriticel  lifting  results  (Fig.  S)  for  e 
NACA  0012  slrfoil  st  M.  -  0.63  and  a  -  2*  indi¬ 
cate  that  the  circulation  model  is  suitable.  The 
same  grids  described  in  the  previous  case  are  again 
used  here.  There  is  a  slight  discrepancy  between 
the  two  codes  in  the  amount  of  leading-edge  suction. 
This  is  most  likely  a  result  of  differences  in  the 
grid,  in  particular,  the  amount  of  grid-point  clus¬ 
tering  around  the  leading  edge.  The  lift  coeffi¬ 
cients  for  TUNA  and  TAIR  were  0.338  and  0.334, 
respectively . 

As  a  fi.ial  two-dimensional  steady-test  case, 
supercritical  lifting-flow  solutions  are  compared 
in  Fig.  6  for  a  NACA  0012  at  H,  ■  0.75  and  0-2°. 

A  finer  grid  with  dimensions  of  151x31  was  used  for 
the  TAIR  computation.  The  two  solutions  are  in 
good  agreement  and  the  shock  profile  for  the  present 
method  la  quite  sharp  despite  the  coarser  grid.  The 
lift  coefficients  for  TUNA  and  TAIR  for  these  two 
computations  were  0.5963  and  0.5881,  respectively. 

The  TAIR  solution  for  each  of  the  airfoil  cases  does 
not  compute  a  stagnation  point  at  the  trailing  edge, 
as  in  the  present  method,  because  the  density  is 
extrapolated  to  avoid  potential  problems  from  the 
grid  singularity.  The  present  method  does  not 
extrapolate  the  density  for  these  computations  but 
It  has  been  found  that  in  some  fine-grid  cases  this 
causes  oscillations  and  overshoots  near  the  trailing 
edge. 

It  is  worth  noting  that  for  supercritical  cases, 
the  AF2  algorithm  has  been  shown  to  provide  faster 
convergence  rates  than  ADI. 24  But  as  previously  ( 
discussed,  the  judicious  addition  of  certain  tem¬ 
poral  damping  terms  may  provide  convergence  rates 
for  ADI  that  approach  AF2.  These  modifications,  as 
yet,  have  not  been  tested.  Convergence-rate  com¬ 
parisons  for  supercritical  flows  with  large  regions 
of  supersonic  flow  and  strong  shock  waves  have 
indicated  that  AF2  converges  2  to  3  times  faster 
than  the  present  ADI  scheme. 

Unsteady  flow  results  are  shown  in  Fig.  7  for 
an  NACA  64A010  airfoil  sinusoidally  oscillating  in 
plunge  +1  about  o  -  0° .  The  reduced  frequency, 
k  -  wc/U„,  is  0.4  and  -  0.80.  The  lift  coeffi¬ 
cient  is  plotted  va  time  during  the  fourth  cycle  of 
oscillation.  Also  shown  are  computed  results  from 
an  Euler  equation  solution.18  The  magnitudes  of 
the  lift  are  in  excellent  agreement,  with  only  about 
a  12  discrepancy.  In  addition,  the  computed  phases 
are  in  good  agreement.  The  Euler  solution  displays 
a  phase  lag  of  approximately  35°  In  comparison  to 
38”  as  predicted  by  TUNA.  Also  shown  are  pressure 
coefficient  distributions  at  three  times  correspond¬ 
ing  to  zero  lift  and  the  maximum  (positive  and 
negative)  values  of  the  lift. 

Three-dimensional,  steady-state  computations 
were  performed  for  a  rectangular  wing  of  aspect 
ratio  6  with  an  NACA  0012  airfoil  section  and  com¬ 
pared  with  the  3D  transonic  computer  code  TWING 
(transonic  wing  analysis)  that  solves  the  steady, 
full-potential  equation  using  the  AF2  scheme.  The 
identical  warped  cylindrical  grid  system  sketched 
in  Fig.  1  was  used  for  the  computations  of  the  com¬ 
puter  codes  TUNA  and  TWING.  The  dimensions  were 
100*10x30  with  5  span  stations  on  the  wing,  which 
is  a  reasonably  coarse  grid  in  the  span  direction. 

A  comparison  of  solutions  for  a  subcrltlcal 
lifting  case  with  M,  •  0.63  and  the  wing  at 
a  •  2*  is  shown  in  Fig.  8.  The  pressure  distribution 


on  the  wing  for  ell  five  span  stations  shows  excel¬ 
lent  agreement  between  the  two  codes.  The  lift 
coefficients  at  the  well  station  for  TUNA  and  TWING 
were  0.250  and  0.253,  respectively,  which  represents 
about  a  25Z  reduction  from  the  two-dimensional 
results  previously  obtained.  TUNA  and  TWING  con¬ 
verged  the  maximuB  residual  two  orders  of  magnitude 
for  this  case  in  85  and  100  Iterations,  respectively. 

The  results  for  a  supercritical  lifting  solu¬ 
tion  with  K.  -  0.75  and  the  wing  at  a  -  2*  ere 
compared  in  Fig.  9.  The  two  solutions  are  in  rea¬ 
sonably  good  agreement.  At  each  span  station,  there 
is  a  slightly  larger  expansion  computed  by  TUNA  from 
the  leading-edge  region  of  the  airfoil  section  to 
the  shock  position,  which  is  located  about  2.52 
farther  upstream  than  that  computed  by  TWING.  This 
discrepancy  may  be  due  to  the  different  ways  in 
which  each  code  determines  the  solution  on  the  wing 
extension.  As  mentioned  before,  TUNA  extrapolates 
and  averages  the  solution  there  while  TWING  solves 
the  full-potential  equation  on  this  section.  The 
shock  profiles  computed  by  TWING  appear  to  be 
sharper  and  the  reexpansion  singularity  is  not 
captured  by  TUNA. 

The  two-dimensional  grids  employed  here  are 
0-type  grids.  Even  though  these  grids  require 
special  attention  at  the  trailing  edge,  a  result  of 
the  mapping  singularity,  they  are  preferred  over 
C-  and  H-type  meshes  because  they  are  more  grid- 
point  efficient,  especially  in  the  far-fleld.  Simi¬ 
larly,  warped  spherical  grids  offer  the  same  advan¬ 
tage  in  grid-point  efficiency  in  three  dimensions 
over  current  wing-type  grids  that  are  typically 
constructed  in  a  manner  similar  to  warped  cylindri¬ 
cal  grids  used  for  the  previous  three-dimensional 
calculations.  The  use  of  spherical  grids  requires 
special  treatment  to  handle  the  mapping  singularity 
on  the  axes.  This  capability  has  been  incorporated 
in  the  present  method. 

To  demonstrate  the  ability  to  employ  spherical 
grids,  the  three-dimensional  incompressible  flow 
over  a  sphere  is  presented  as  a  simple  check  case. 

A  comparison  between  the  exact  incompressible  flow 
solution  and  the  present  method  with  M„  »  0.01  is 
shown  in  Fig.  10.  A  40*21*22  radially  stretched 
spherical  grid  was  used.  The  solution  shown  is  the 
pressure  distribution  on  the  upper  surface  of  the 
sphere  from  leading  to  trailing  edge  in  the  plane  of 
symmetry  perpendicular  to  the  axes.  The  agreement 
is  correct  to  plottable  accuracy.  The  maximum 
residual  converged  six  orders  of  magnitude  in 
30  iterations  for  this  case. 

A  more  complicated  geometry  using  a  spherical 
grid  topology  is  an  ellipsoid  type  wing  with  a  chord 
of  1,  aspect  ratio  of  4,  and  thickness  of  1/4. 

Here,  S6  with  the  sphere,  the  axes  of  the  spherical 
grid  are  aligned  with  the  y  axis  and  the  flow  is 
in  the  x-direction.  The  surface  grid-point  distri¬ 
bution  and  other  partial  views  of  the  grid  are  shown 
in  Fig.  11.  .he  pressure  coefficient  on  the  upper 
surface  in  the  y  -  0  plane  is  shown  in  Fig.  12  for 
M„  ■  0.77.  A  relatively  strong  shock  can  be  seen  at 
about  852  chord.  The  shock  profile  is  not  particu¬ 
larly  sharp  but  the  grid  is  fairly  coarse  with  only 
60  points  around  the  body  in  the  1  direction.  The 
residual  converged  3  orders  of  magnitude  in 
100  iterations  for  this  case. 

The  current  version  of  TUNA  has  been  programmed 
and  run  on  the  CRAY  IS.  The  program  requires 
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approximately  0.0001  sac  of  CPU  time  par  grid  point 
par  itaration.  The  coda  la  presently  in  an  Ineffi¬ 
cient  state  and  no  effort  has  bean  made  to  taka 
advantage  of  the  vector  capabilities  of  tba  CRAY . 


V.  Concluding  Remarks 

An  Implicit,  approximate-factorisation,  finlte- 
dlffaranca  algorithm  has  been  developed  for  solving 
the  conservative  full-potential  equation  for  com¬ 
puting  steady  and  unsteady  transonic  flows  In  two- 
and  threa-dimenalone .  A  general  coordinate  trans¬ 
formation  was  employed  for  the  simple  and  accurate 
treatment  of  arbitrary  body  surfaces.  This  algo¬ 
rithm  has  been  programmed  Into  a  computer  code 
named  TUNA.  Results  have  been  presented  for  several 
two-  and  three-dimensional  steady  and  two- 
dimensional  unsteady  flows  to  verify  the  secure  '.y 
of  the  algorithm. 
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Fig.  1  Schematic  of  cylindrical  grid  system  used  , 
for  wing  calculations  employing  a  plane  of  symmetry. 


Fig.  2  Two-dimensional  pressure  coefficient  distri¬ 
bution  comparison:  NACA  0012  airfoil,  M,  •  0.72, 
a  -  0*. 


a*Bolat,  T.  L. ,  and  Ballhaus,  w.  F.,  tu., 
Conservative  8chcmaa  for  the  Full  Potential  Equa¬ 
tion  Applied  to  Transonic  Flows,"  AIAA  Journal. 
Vol.  17,  Fab.  1979,  pp.  145-152. 


Fig.  3  O-mesh  type  finite  difference  grid  for  an 
NACA  0012  airfoil  numerically  generated  using  GRAPE 
code  (100*31  grid  points). 


NUMBER  OF  ITERATIONS 

Fig.  4  Residual  history  comparison  between  the  time 
accurate  algorithm  and  ADI  relaxation  option: 

NACA  0012  airfoil,  N.  -  0.72,  a  -  0*. 
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Fig.  6  Two-dtnenatona 1  pressure  coefficient  distri¬ 
bution  coaparlson:  NACA  0012  airfoil,  M„  ”  0  71, 
a  -  2*. 


Fig.  7  Time  variation  of  the  lift  coefficient  and  various  pressure  coefficient  distributions  for  an 
NACA  64A010  airfoil  sinusoidally  oscillating  in  plunge,  K„“  0.80,  k  -  0.40,  Op  •  1*  sin(kH»t). 


1.0 


Fig.  10  Comparison  of  computed  pressure  coefficient  with  Incompressible  flow  theory  for  flow  about  a 
about  a  sphere. 


Fig.  11-  Spherical  grid  used  for  ellipsoidal  wing,  a)  Flanform  view  of  surface  grid  point  distribution; 
Partial  grid  views  of  b)  Symmetry  plane  (z  •  0) ,  c)  Midspan  symmetry  plane  (y  *  0) ,  and  d)  Midchord 
symmetry  plane. 
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A  ZONAL  APPROACH  FOR  THE  STEADY  TRANSONIC 
SIMULATION  OF  INVISCtD  ROTATIONAL  FLOW 


Neal  M.  Cbaderjiaor  and  Joseph  L.  Steger* 
Stanford  Univtrtiiy,  Stanford,  California 


Abstract 

A  finite  difference  tonal  method  is  developed  to 
compute  steady  inviscid  transonic  flow  by  coupling  a 
semi-flux  split  form  of  the  Euler  equations  in  a  vor- 
ticity  producing  zone  with  a  tone  of  scalar  and  vec¬ 
tor  (i.e.,  dual)  potential  equations.  The  dual  potential 
equations  permit  vorticity  convection,  but  not  produc¬ 
tion,  and  are  efficiently  solved  as  an  iteratively  decou¬ 
pled  set  of  scalar  equations.  Zonal  results  presented  for 
a  nonlifting  biconvex  airfoil  on  a  stretched  Cartesian 
grid  show  substantial  savings  in  CPU  time  compared 
to  solving  the  semi-flux  split  Euler  equations  alone. 
The  dual  potential  equations  also  provide  an  alternate 
way  of  treating  potential  flows  with  circulation.  This 
has  been  demonstrated  by  computing  a  subcritical  flow 
over  a  lifting  airfoil  using  generalized  curvilinear  coor¬ 
dinates. 


L  Introduction 

The  development  of  efficient  numerical  solution 
techniques  for  simulating  rotational  compressible  flow 
has  always  been  subject  to  two  competing  philosophies. 
In  one  approach  the  Navier-Stokes  or  Euler  equations 
are  programmed  throughout  the  entire  domain.  The 
overall  computer  code  is  quite  general  in  its  applicability 
and  the  programming  is  straight  forward  in  the  sense 
that  only  one  equation  set  is  coded.  In  the  other  ap¬ 
proach  the  flow  field  is  partitioned  into  zones  which  use 
the  simplest  and  most  efficiently  solved  equations  pos¬ 
sible.  For  example  the  Navier-Stokes  equations  m3y 
be  used  in  one  zone,  the  nonlinear  potential  equations 
in  another,  and  perhaps  linear  theory  in  the  remain¬ 
ing  field.  This  zonal  approach  has  the  potential  for 
saving  substantial  amounts  of  the  computer  resource. 
However,  a  zonal  code  is  significantly  more  complex 
to  program  because  several  computer  codes  must  es¬ 
sentially  be  written  (one  for  each  governing  equation 
set)  and  each  zone  must  be  interfaced.  Moreover,  con¬ 
siderable  care  must  be  taken  in  interfacing  zones  be¬ 
cause  poor  interfacing  can  result  in  computational  in¬ 
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efficiency  thus  offsetting  the  advantages  of  a  zonal  ap¬ 
proach. 

Whether  the  use  of  a  zonal  method  will  ultimately 
prove  better  than  just  using  a  single  general  equation 
set  Is  a  question  that  will  likely  remain  moot,  but 
several  developments  have  occurred  which  In  our  view 
make  the  zonal  approach  more  attractive  than  pre¬ 
viously.  The  first  development  Is  that  finite  difference 
lnvlsctd  and  viscous  boundary  layer  Interaction  algo¬ 
rithms  are  becoming  available  for  solving  separated 
viscous  flow  and  these  schemes  appear  to  be  much 
faster  than  Navier-Stokes  algorithms.  The  second  de¬ 
velopment  Is  the  appearance  of  efficient  methods  for  In¬ 
cluding  rotatlonallty  effects  Into  essentially  potentlal- 
ilke  governing  equation  sets.  It  Is  this  second  develop¬ 
ment  which  Is  of  Interest  to  us  In  this  paper. 

A  zonal  algorithm  becomes  more  advantageous 
If  potential-Uke  equations  can  be  generalized  to  simu¬ 
late  rotational,  nonlsentroplc  flow  without  loss  of  com¬ 
putational  efficiency.  This  Is  because  more  of  the  flow 
fleid  can  be  treated  with  the  more  efficient  generalized 
potential  code  and  fewer  zones  and  thus  fewer  Interface 
boundaries  have  to  be  Introduced. 

In  this  paper  we  study  a  two  zone  method  for 
solving  Inviscid  transonic  flow  without  requiring  the 
flow  to  be  Irrotatlonal.  Shock  waves  are  captured  In 
the  first  zone  using  a  semi-flux  split  implicit  finite 
difference  method  to  solve  the  Euler  equations  in  strong- 
conservation-law  form.  In  this  algorithm  flux  splitting 
Is  used  In  the  direction  along  the  body,  and  central 
differencing  Is  used  In  the  direction  away  from  the 
tody.  Because  of  this  differencing  structure,  this  al¬ 
gorithm  Is  readily  extended  to  Incorporate  thin  layer 
viscous  terms  although  we  do  not  exercise  this  op¬ 
tion  here.  A  dual  potential  formulation  Is  solved  In 
a  zone  away  from  the  shock  and  Is  able  to  correctly 
convert  but  not  generate  vorticity.  This  second  zone 
slightly  overlaps  the  shock  zone  and  extends  over  the 
remainder  of  the  flow  field.  The  dual  potential  scheme, 
w  hich  Is  similar  to  a  formulation  due  to  Hafez  and 
Lovell1,  decomposes  a  velocity  fleid  Into  gradients  of 
scalar  potential  and  vector  potential.  The  vector  poten- 
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tlal  function  accounts  for  vorticlty.  and  as  a  conse¬ 
quence.  docs  not  require  a  circulation  cut  for  a  lifting 
airfoil  as  does  a  scalar  potential  formulation.  The  cru¬ 
cial  advantage  of  the  dual  potential  formulation  over 
a  primitive  variable  formulation  for  the  Euler  equa¬ 
tions  Is  that  the  resulting  dual  potential  equations  are 
weakly  coupled.  They  can  thus  be  treated  In  scalar 
mode  and  can  take  advantage  of  all  of  the  numerical 
efficiencies  developed  for  the  transonic  potential  equa¬ 
tion. 

In  the  following  sections  we  discuss  the  govern¬ 
ing  equation  sets  and  zonal  coupling  concepts.  A  set  of 
boundary  conditions  are  developed  which  permit  vor- 
tldty  convection,  and  a  new  way  of  treating  circula¬ 
tion  by  taking  advantage  of  the  stream  function-like 
properties  of  the  dual  potential  equations  are  presented. 
Finally,  the  numerical  algorithms  are  developed  fol¬ 
lowed  by  a  discussion  of  results  and  concluding  remarks. 

D.  Zonal  Formulation 
Governing  Equations 

The  Sow  Seld  has  been  zoned  between  two  equa¬ 
tion  sets.  In  regions  of  expected  vorticlty  generation, 
here  shock  waves,  the  Euler  equations  are  solved  In- 
strong-conservation-law  form.  In  transformed  coor¬ 
dinates  these  equations  are  given  by 

+  dqCi  ■>  0  (1) 

where 

G  "  +  «!»<*)/  J 

and 

('£')  •  *-(/;,) 

V.u(e  +  p)/  Vv(e  +  py 

In  the  above  equations  p  Is  the  density,  u  and  v  the 
Cartesian  velocities  In  the  x  and  y-dlrectlons.  p  = 
<7—  1  )|e  —  Is  the  pressure  and  e  the  total  energy 
per  unit  volume.  The  fluid  speed  Is  q.  The  Cartesian 
flux  vectors  F  and  G  correspond  to  the  r  and  in¬ 
directions.  and  J  Is  the  determinant  of  the  transfor¬ 
mation  Jacobian. 


The  variables  have  been  nondlmenslonallzed  by 
the  free  stream  density  and  velocity  as 

P  **  P/P oo 
1  —  a/«oo 
*  —  »/«<» 

P  —  P/(Poo«2o) 

*  —  C/(POO«*50) 


where  the  —  has  been  suppressed  for  convenience. 

la  the  preseat  test  program  we  use  small  distur¬ 
bance  thin  airfoil  boundary  conditions  so  only  simple 
stretching  transforms  are  needed  to  duster  grid  points. 
The  stretchings  are  of  the  form 

<  -  «*) 

*1  —  u(y)  (2) 


Consequently, 

t,  —  o  .  if,  —  0 

and  Eq.  (1)  Is  very  much  like  Its  Cartesian  counterpart. 

In  the  remaining  part  of  the  flow  fleld  a  restricted 
form  of  the  steady  Euler  equations  are  used.  In  non- 
dimensional  variables  these  are  given  by 
continuity 

fi,(pu)  +  dy(pv)  «  0  (3) 

Crocco  (vorticlty)  equation 
9,v  —  dvu  «=»  -(7Jlda)”,(»d,e  —  udyt)  (e) 
constant  entropy  along  a  streamline 

Ud,«  +  VdyS  “  0  (3) 

Bernoulli  equation 

P^fl  +  ^—Ad^t-tt’-v*)]  '  e—  (61 


The  nondlmenslonallzed  variables  are  consistent  with 
those  used  for  the  conservation-law-form  equations 


with  nondiraensional  entropy  i  m,  (a  —  s<x>)/R.  Once 
again  the  has  been  suppressed  with  R  the  perfect 
gas  constant,  7  the  ratio  of  specific  heats,  and  M  the 
Mach  number. 

In  regions  of  entropy  production,  for  example 
across  a  shock  wave  or  in  viscous  layers,  Eq.  (5)  is 
invalid.  The  simplified  equations  as  given  also  assume 
uniform  stagnation  enthalpy,  although  it  is  relatively 
easy  to  incorporate  this  effect  into  the  equations  and 
to  solve  as  well  the  equation 

ud,h,t  +  vduh,t  —  0  (7) 

In  order  to  take  advantage  of  efficient  poten¬ 
tial  flow  solvers,  a  dual  potential  representation  of 
velocity  is  introduced  into  the  simplified  equation  set. 
Specifically,  the  velocity  components  are  decomposed 
into  scalar  potential  and  vector  potential  functions  which 
for  two  dimensions  are  given  by 

u  9,4  +  9vt/>  *=  da  + 

V=»  9v4~-d*4  =*  da  (8)  ' 


where  v  is  the  third  vector  potential  component.  Intro¬ 
ducing  the  dual  potential  velocity  relations  into  Eqs. 
(3)  and  (4)  we  obtain' 

continuity 

djMda+*»)l  +  9wMd,-da)l  —  0  (9) 

Crocco  or  vorticity 

•*  (»«,  -  u«v)  ■»  -n  (10) 

And  finally,  with  transformation  into  general  coor¬ 
dinates 


d((patd(  +  pojd^  I  +  #n(pa2d<  +  t><*s4n) 
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where 
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and 

5<(a,d<  +  oade)  +  *e(oratf<  +  asd,) 

—  (7J +  *,*') +  ?***)! 

(12) 

where 

«  “  0*e  +  n,49  +  01>t  +  ntf9 

v  ■«  (,dt  +  q»de  -  Ol>t  -  *U4>9 


For  the  current  zona)  algorithm  applications  we 
will  again  restrict  the  transforms  to  the  simple  stretch¬ 
ings  £  —  £(x)  and  n  mm  jj(y).  However,  the  fully  trans¬ 
formed  dual  potential  equations  alone  will  be  used  to 
solve  a  subcritical  lifting  airfoil  flow  in  order  to  il¬ 
lustrate  their  applicability  to  flows  with  circulation 
without  the  need  to  impose  special  cuts  in  the  field. 
Zone  Partitioning 

As  sketched  in  Fig.  la  and  lb,  the  conservation- 
law-form  equations  of  mass,  momentum,  and  energy 
are  to  be  solved  in  shock  regions  where  vorticity  is 
generated.  Ideally  this  conservation-law  zone  would 
be  kept  as  small  as  possible  as  shown  in  Fig.  la.  This, 
however,  will  require  shock  wave  tracking  logic  and 
consequently  we  are  currently  defining  a  much  more 
generous  stationary  zone  for  Eq.  (1)  as  sketched  in 
Fig.  lb.  In  the  remaining  flow  field  the  simplified 
governing  equations  are  solved  in  terms  of  the  dual 
potential  variables.  For  uniform  incoming  flow  entropy 
correction  terms  are  only  needed  behind  the  shock 
wave. 

Boundary  Conditions 

The  inviscid-conservation-law  equations  have  been 
restricted  to  a  zone  about  the  expected  shock  wave 
which  overlaps  a  larger  zone  governed  by  the  dual 
potential  equations  (see  Fig.  lb).  Except  along  the 
body  surface,  all  boundary  conditions  for  the  conserva¬ 
tion-law  zone  can  be  supplied  from  the  overlapping 
dual  potential  zone.  Along  the  body  surface  we  impose 
the  tangency  condition 

v  "3) 
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and  use  linear  extrapolation  to  supply  the  three  remain¬ 
ing  variables.  A  thin  airfoil  approximation  is  used  by 


imposing  these  boundary  conditions  at  y  »  0  rather 
than  on  the  body  surface.  This  boundary  treatment 
permits  verification  of  the  method  without  unduety 
complicating  the  computer  code.  This  will  also  lead  to 
some  solutioa  discrepancies  when  results  are  compared 
to  more  exact  theories. 

The  dual  potential  equations  resolve  the  flow  in 
the  remaining  domain.  Representing  the  velocity  by 
derivatives  of  scalar  functions  allows  some  freedom  in 
the  choice  of  boundary  conditions.  We  will  therefore 
describe  the  dual  potential  boundary  conditions  for 
our  specific  application  which  is  the  transonic  flow  over 
a  symmetric  airfoil  at  zero  angle  attack  (see  Fig.  2). 

In  prescribing  a  set  of  boundary  conditions  for 
the  potential  functions,  we  adopt  the  point  of  view  that 
0  js  a  perturbation  on  0  in  the  sense  that  0  represents 
an  irrotational  flow  and  0  will  only  be  non-zero  if  there 
is  vorticity  or  lift.  For  uniform  incoming  flow  about 
a  symmetric  thin  airfoil  an  appropriate  set  of  far  field 
boundary  conditions  for  0  are  given  by 

*-*  (14) 

The  tangency  condition  is  imposed  at  y  «=  0  by 


V»*  +  (0*  + 
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An  appropriate  upstream  far  field  boundary  condition 
for  v  is  given  by 


Note  that  the  0-function  is  zero  at  the  upper  left 
corner,  but  can  be  non-zero  at  the  upper  right  corner 
provided  there  is  vorticity  in  the  airfoil  wake.  This 
problem  is  resolved  by  a  limiting  process  on  the  upper 
boundary  where 


,r  t.d* 

*  — OO 


constant 


{£«*«- 0  (2») 

As  the  upper  boundary  is  extended  farther  away  from 
the  airfoil  in  both  the  x  and  y-directions,  the  0-function 
variation  will  be  small  yet  permit  a  change  from  zero 
at  the  upper  left  boundary  to  the  appropriate  value 
on  the  upper  right  boundary.  This  process  has  been 
numerically  verified  by  observing  that  the  top  bound¬ 
ary  velocities  approach  the  uniform  condition. 

Equation  (5)  is  a  convective  entropy  equation  that 
can  be  marched  in  the  x -direction.  Initial  entropy 
data  behind  the  shock  is  obtained  from  the  Euler  equa¬ 
tions.  The  entropy  is  constant  on  the  airfoil  and  sym¬ 
metry  streamline  in  the  wake  region. 

Zone  Interfacing 

The  two  different  sets  of  governing  equations  are 
interfaced  at  their  zone  boundaries.  In  Fig.  3,  for  ex¬ 
ample,  boundary  values  can  be  applied  to  the  conserva¬ 
tion  equations  (1)  along  the  curve  efgh  by  differentiating 
0  and  0.  Along  the  inner  boundary  curve  abed  one  can 
specify  «,  u,  and  v  from  the  solution  of  Eq.  (1)  so  that 


The  ^-function  behaves  much  like  a  stream  function 
which  motivates  the  boundary  condition  on  the  lower 
boundary  (y  «=  0) 


0»  +  0»  "  u»r«ci/(«4 
0|(  ~  ™  VtfMCf'/feS 


If  this  were  a  lifting  airfoil  problem  then  0  would  be 
chosen  as  some  non-zero  constant  on  the  airfoil  surface 
in  order  to  satisfy  a  Kutta  condition  (see  Section  III). 
A  shock  induced  rotational  flow  can  give  rise  to  a 
velocity  defect  in  u  at  the  right  boundary  so  that  0B 
must  be  free  to  vary.  A  suitable  downstream  far  field 
boundary  condition  to  ensure  that  v  0  for  0  =  r  is 
given  by 


On  the  top  boundary  we  set  u 
x,  we  impose 


I,  and  because  0  «= 


supply  derivative  boundary  conditions  for  0  and  0.  By 
overlapping  the  boundaries  as  illustrated,  information 
can  be  efficiently  transferred  from  one  domain  to  the 
other.  In  overlap  regions  both  equation  sets  are  being 
solved  and  are  differentially  equivalent  if  the  shock 
wave  is  avoided 

We  have  used  another  way  to  interface  the  two 
zones  together.  From  the  vantage  point  of  computa¬ 
tional  simplicity,  it  is  desirable  to  interface  the  equa¬ 
tions  using  as  little  special  logic  as  possible.  Except 
across  the  shock  wave,  the  dual  potential  equations  are 
everywhere  equivalent  to  Eq.  (I)  for  steady  flow  with 
a  uniform  incoming  stream.  Even  across  the  shock 


Eq.  (9)  is  valid  in  the  sense  that  it  perservea  the  cor¬ 
rect  weak  solution.  Consequently  we  can  solve  Eq. 

(9)  throughout  the  entire  domain  without  having  to 
define  a  hole  with  special  boundary  condition  treat¬ 
ment  as  discussed  above.  This  simplifies  the  computer 
code  logic,  especially  since  we  avoid  solving  Eq.  (22) 
on  the  zonal  boundary.  Also,  by  avoiding  the  hole, 
we  solve  the  continuity  prediction  equation  for  d  more 
implicitly.  That  is,  we  avoid  iteratively  lagged  inner 
boundary  data  which  must  be  updated  from  the  solu¬ 
tion  of  the  conservation-law  form  equations. 

The  Crocco  equation  which  is  used  as  a  predic¬ 
tion  equation  for  d  does  not  admit  a  Rankine-Hugoniot 
jump  solution.  Nevertheless,  we  have  also  used  this 
equation  throughout  the  entire  domain,  including  across 
the  shock  so  as  to  avoid  integrating  Eq.  (22)  at  the 
zone  boundary  and  to  improve  implicitness.  However, 
we  supply  the  right-band-side  vorticity  function  of  Eq. 

( 10)  in  the  zone  abed  directly  from  the  solution  of  Eq. 
(1).  That  is,  we  solve  (shown  in  Cartesian  form) 

■  (Uy  —  Vj)(f nei/Ud  (23a) 

or 

d»«  +  *■  (tA/*)— '(«*»  -  ueJ.Me< (236) 

where  s  or  derivatives  of  u  and  v  are  specified  from  the 
conservation-law  form  equations.  While  this  equation 
is  not  strictly  valid  across  the  shock,  we  have  used  it 
together  with  Eq.  (10)  so  as  to  avoid  coding  an  inner 
hole  boundary  condition  (i.e.,  on  abed  in  Fig.  3)  and 
to  improve  implicitness. 

The  equation  for  convection  of  entropy,  Eq.  (5), 
cannot  be  used  across  the  shock.  Consequently,  in  the 
entropy  production  zone,  here  bounded  by  the  curve 
abeda  of  Fig.  3,  we  update  entropy  from  the  solution 
of  Eq.  (1)  using  the  thermodynamic  relation 

«"(*»-  i)"1  <"b(7  - 1) Af2o(e  “  £P?2)M  (24) 

Summarizing  these  concepts,  we  have  opted  for 
the  following  interface  scheme.  For  the  vorticity  produc¬ 
tion  zone  enclosed  by  the  curve  efghe  shown  in  Fig.  3, 
Eq.  (1)  is  solved.  Values  of  p,  pa,  pv,  and  e  are  supplied 
on  the  outer  boundary  of  this  zone  from  the  solution 
of  the  dual  potential  equations.  The  dual  potential 
equations,  Eqs.  (9)  and  (10),  are  solved  over  the  entire 
flow  field  domain.  However,  within  the  vorticity  zone 
circumscribed  by  abeda,  Eq.  (23)  is  used  in  place  of 
Eq.  (10).  This  equation  substitution  or  modification 


is  readily  accomplished  with  simple  flags  and  only  one 
bask  set  of  computer  algorithms  are  needed  to  solve 
the  dual  potential  equations  over  the  entire  domain. 
In  this  same  zone  Eq.  (24)  replaces  Eq.  (5). 

HI.  Circulation  Treatment 

In  a  code  based  solely  on  a  scalar  potential  it  is 
necessary  to  build  cuts  into  the  field  in  order  to  treat 
flow  with  circulation.  Across  these  cuts  the  scalar  (i.e., 
velocity)  potential  is  discontinuous  and  the  difference 
equations  must  be  specially  coded  to  account  for  this 
jump.  Because  the  dual  potential  equations  allow  vor¬ 
ticity  to  convect  through  the  flow  field,  it  is  no  longer 
necessary  to  build  such  circulation  cuts  into  the  field. 
This  simplicity  partially  compensates  for  having  to 
solve  Eq.  (10). 

In  order  to  treat  lifting  airfoils  with  the  dual 
potential  equations  we  adjust  on  the  airfoil  so  as  to 
satisfy  the  Kutta  condition.  In  essence  we  handle  the 
vector  potential  component  i>  much  as  if  it  were  the 
«tream  function.  When  solving  a  flow  with  a  stream 
function  formulation,  the  value  of  the  stream  function 
on  the  body  is  a  constant  which  is  determined  by  the 
Kutta  condition.  Here  we  implement  the  Kutta  con¬ 
dition  by  requiring  that  the  magnitude  of  the  trailing 
edge  velocities  match  on  the  upper  and  lower  surfaces 
(see  Fig.  4).  That  is 

"*  Qlowor  (25) 

or 


For  isentropie  flows  this  is  equivalent  to  pmpp,r  •“  Plow 
The  Cartesian  velocity  components  u  and  v  can  be 
eliminated  in  terms  of  and  as  given  by  Eq.  (12). 
Making  use  of  the  boundary  condition  that  i<  is  con¬ 
stant  on  the  body  (rpf  *■  0),  Eq.  (26)  becomes 


(27) 

This  relation  can  be  differenced  to  provide  the  value 
of  $  on  the  body  that  satisfies  the  Kutta  condition. 


IV.  Numerical  Algorithms 
Implicit,  approximately  factored  (AF)  Aoite  differ¬ 
ence  schemes  are  used  to  solve  the  governing  equations. 

A  semi-flux  spilt  algorithm  Is  used  for  Eq.  (1).  and  al¬ 
ternating  direction  methods  are  used  with  Eqs.  (11) 
and  (12). 

Semi-Flux  Split  Algorithm 

The  semi-flux  spilt  Euler  equations  In  strong  eon- 
servatlon-law-form.  Eq.  (1).  are  solved  In  delta  form 
using  the  Implicit  approximately  factored  flnlte  difference 
scheme 


|/  +  AV«((,A+)  +  X 

(/  +  -  —hR*  (28) 


where 

r  -  (5*(^J/“,r+)  +  sfaj-'F-HtihiJ-'Gi  | 
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and  F*  are  flux  spilt  vectors  as  derived  by  Steger  and 
Warming*  while  A*  are  the  Jacobian  matrldes  °fi£. 
The  vector  of  primitive  variables  Is 


The  parameter  A  Is  a  constant  pseudo-time  step 
chosen  to  accelerate  the  convergence  rate.  All  metric 
quantities  are  computed  with  second-order  accurate 
difference  formulas  using  one  sided  differences  at  the 
flow  boundaries  and  central  differencing  Interior  to  the 
Row  boundaries.  The  difference  operators  In  Eq.  (28) 
are  defined  by 
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where  A(  »  A»j  ■■  I.  For  simplicity,  Qrst-order 
accurate  one  sided  derivative  operators  are  used  In 
the  Implicit  part  and  second-order  accurate  one  sided 
differences  are  used  In  the  residual.  The  converged 
solution  will  be  second-order  accurate  and  uncondi¬ 
tional  linear  stability  Is  still  retained.  The  flux  split¬ 
ting  In  the  (-direction  Is  Inherently  dissipative:  however 
fourth  order  dissipation  Is  added  In  the  q-dlrecilon 
with  t~0(A). 

The  first  factor  ofEq.  (28)  Is  solved  by  sweeping 
forward  In  (  and  Inverting  block  trldlagonal  systems  of 
equations  In  the  q-dlrectlon.  The  second  factor  forms 
an  upper  bldlagonal  matrix  which  Is  solved  by  a  simple 
back  sweep.  In  supersonic  regions  the  second  factor 
disappears  (F~  =  0)  so  the  AF  scheme  reduces  to  a 
single  direct  Inversion  of  the  ttme-llnearlzed  equations. 
In  our  application  the  flow  In  the  vortldty  produc¬ 
tion  zone  Is  mostly  supersonic;  consequently,  the  AF 
scheme  Is  very  efficient  In  that  zone. 

Dual  Potential  Algorithms 

The  dual  potential  equations  are  weakly  coupled 
In  the  sense  that  one  can  effectively  solve  Eq.  (12)  for 
0  as  a  function  of  lower  derivative  terms  that  are  fixed 
at  a  previous  Iteration  level.  Once  Is  predicted.  Eq. 
(11)  serves  as  a  prediction  equation  for  4>.  At  this  point 
entropy  can  be  obtained  from  Eq.  (8)  and  density  Is 
updated  from  the  Bernoulli  equation. 

Implicit  approximately  factored  finite  difference 
algorithms  are  used  to  solve  the  dual  potential  equa¬ 
tions.  We  first  describe  the  numerical  algorithm  for 
Eq.  (11)  which  Is  very  similar  to  what  Is  used  for 
the  transonic  full  potential  equations.  Our  Implemen¬ 
tation  follows  that  of  Sieger  and  Caradonna4.  The 
AF  scheme  Is  second-order  accurate  In  subsonic  flow 
regions  and  Is  first  or  second-order  accurate  In  super¬ 
sonic  flow  regions.  Upwlndlng  In  the  supersonic  region 
Is  accomplished  by  using  Holst's  upwind  shifted  den¬ 
sity  scheme*. 

The  AF  scheme  for  Eq.  (11)  Is  given  In  delta 
form  as 


! I  —  AV,(p"n*  J~l  X 
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where 
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Tbe  j,  k-indkes  correspond  to  the  {,  17 -directions  and 
are  suppressed  except  to  indicate  midpoint  values.  The 
pseudo-time  step  A  and  tbe  relaxation  parameter  u  are 
chosen  so  as  to  accelerate  the  iterative  convergence. 
The  {-difference  operators  are  defined  by 
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Equation  (29)  is  solved  in  the  usual  fashion  by  invert* 
ing  scalar  tridiagonal  systems  of  equations  in  first  the 
9  and  then  the  {-directions. 

A  second-order  accurate  difference  for  the  tan- 
gency  boundary  condition,  Eq.  (IS),  is  implemented 
using  central  differencing  for  z-derivatives  and  three 
point  one  sided  differencing  for  tbe  y -derivatives.  A 
tridiagonal  solution  provides  d  at  the  lower  boundary 
with  the  terms  evaluated  explkitly. 

A  similar  algorithm  for  the  vorticity  equation  is 
used  to  update  i>. 


where  again  for  convenience  A{  ™  I.  Similar  expres* 
sions  define  the  q-difference  operators.  Stability  in  su¬ 
personic  regions  is  maintained  by  shifting  the  density 
upwind  in  the  {-erection  according  to  the  formula 
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with 


v,  -  max[0,(M$  -  1)C)  ,  1<C<2 


At  subsonic  points,  t/j  »  0  and  Vj+i/t  is  second-order 
accurate.  The  parameter  9  »  2  provides  second-order 
accuracy  at  supersonk  points  and  9  »  0  gives  first- 
order  accuracy.  In  practice  a  value  of  9  ■■  1.8  is  used 
for  Bows  with  weak  shocks  and  9  is  dropped  to  sero  for 
stronger  shock  flows.  The  term  ffV((Ad”)  provides 
some  additional  iterative  damping  and  compensates 
somewhat  for  not  properly  linearizing  the  density  term 
for  inclusion  on  the  left  hand  side  of  Eq.  (29).  Tbe 
metrics  are  computed  by  the  following  relations  in 
order  to  capture  uniform  flow 
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Here  though  the  equation  remains  elliptic  throughout. 
Tbe  Neumann  boundary  conditions  Eqs.  (18)  and  (19) 
provide  \l>  on  the  downstream  and  top  boundaries  using 
second-order  accurate  one  sided  differences. 

Finally,  Eq.  (5)  is  solved  for  entropy.  Tbe  stretched 
Cartesian  grid  used  for  the  zonal  test  problem  is  closely 
aligned  with  the  flow,  so  it  is  a  simple  matter  to  ob¬ 
tain  an  update  of  t  by  marching  the  equation  in  the 
{-direction.  A  finite  difference  scheme  for  marching 
in  {  is  given  by 


f\«  + 
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where  is  a  three  point  backward  difference  in  (  and 
j,  is  a  central  difference  in  q.  At  each  {-station  a 
scalar  tridiagonai  inversion  in  q  is  required. 

For  lifting  airfoils  we  must  also  satisfy  the  Kutta 
condition  by  differencing  Eq.  (27)  and  solving  for  0 
on  the  body.  Tbe  {-derivatives  are  central  differenced 
and  tbe  q-derivstives  are  forward  differenced.  Solving 
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— •  second  order 


The  j.k-indices  correspond  to  the  O-grid  topology  in 
Fig.  4.  The  Kutta  condition  is  evaluated  explicitly  at 
the  end  of  each  iteration  of  Eqs.  (29)  and  (30). 

V.  Results  end  Discussion 

A  nonlifting  biconvex  airfoil  was  chosen  as  a  test 
problem  to  verify  the  zonal  approach  as  an  alternative 
to  the  Euler  equations  for  solving  transonic  rotational 
flow.  The  subcritical  flow  about  a  NACA  0012  air* 
foil  at  angle  of  attack  was  used  to  demonstrate  the 
capability  of  the  dual  potential  equations  to  treat  lift 
without  the  use  of  circulation  cuts.  All  computations 
were  performed  on  a  C DC  7600. 

The  biconvex  airfoil  problem  was  computed  on 
a  stretched  Cartesian  grid  with  the  small  disturbance 
boundry  condition  Eq.  (15)  imposed  at  the  lower  bound* 
ary,  y~0.  Figure  5  shows  the  grid  in  the  vicinity  of  the 
airfoil  with  the  far  field  boundary  3  chord  lengths  away 
to  either  side  of  the  airfoil  and  5  chord  lengths  away  in 
the  y-direction.  Exponential  stretching  was  used  away 
from  the  airfoil  in  both  the  x  and  y*directions  while 
a  spline  function  distributes  points  along  the  airfoil 
(0<x<l).  This  allows  grid  clustering  at  the  leading 
edge,  trailing  edge  and  the  shock. 

Results  for  a  nonlifting  10%  thick  biconvex  air¬ 
foil  at  Afoo  **  85  are  presented  for  the  flux  split, 
zonal  and  potential  schemes.  Figure  6  compares  the 
flux  split  Cp  solution  with  a  central  differenced  Euler 
solution*.  A  monotone  solution  at  the  shock  for  the 
flux  split  algorithm  is  obtained  by  conservatively  switch¬ 
ing  to  first-order  accuracy  at  a  few  points  bracketing 
the  shock.  This  accounts  for  the  smearing  at  the  foot 
of  the  shock.  A  comparison  of  the  zonal  Cp  solution 


with  the  same  Euler  solution  is  given  in  Fig.  7.  Except 
at  the  shock,  the  flux  split  equations  are  second-order 
accurate  and  require  two  levels  of  zonal  boundary  data 
in  the  {-direction.  The  zonal  steady  state  solution  ac* 
curacy  is  not  sensitive  to  the  amount  of  zonal  overlap. 
Here  we  have  nsed  three  points  of  overlap  in  both  the 
(  and  q -directions.  As  shown  in  Fig.  7,  the  shock 
wave  locations  are  coincident  within  one  grid  point. 
To  ensure  proper  shock  position  it  was  found  to  be 
important  to  zero  out  a  smalt  value  of  entropy  that 
the  flux-split  scheme  produces  before  the  shock;  other¬ 
wise,  the  density  predicted  from  the  Bernoulli  equa* 
tion  is  inaccurate  at  the  shock  and  the  shock  moves 
upstream  two  or  three  grid  points.  Figure  8  shows 
a  comparison  between  the  zonal  and  scalar  potential 
solutions.  Vorticity  effects  are  evident  by  the  zonal 
shock  location  upstream  of  the  irrotational  potential 
shock.  In  these  calculations  we  have  opted  for  first- 
order  accuracy  in  supersonic  regions  by  evaluating  the 
shifted  density  with  0  0.  The  1^  norm  of  the  6- 

'esidual  is  given  in  Fig.  9  for  the  scalar  potential 
solution.  We  remark  that  the  present  computer  code 
has  not  been  optimized  and  uses  constant  pseudo-time 
steps.  The  convergence  history  for  the  zonal  solution 
is  given  in  Fig.  10.  The  convergence  rate  appears  to 
be  limited  by  Eq.  (29).  For  a  zonal  computation  in 
which  the  conservation-law  zone  was  allocated  20% 
of  the  total  number  of  grid  points,  the  zonal  method 
takes  60%  less  CPU  time  per  iteration  than  the  flux 
split  algorithm.  Moreover,  the  zonal  algorithm  only 
required  50%  -  75%  as  many  iterations  that  the  flux 
split  algorithm  does  for  plottable  accuracy. 

The  flow  about  a  12%  thick  biconvex  airfoil  at 
Afoo  ™  -88  generates  a  solution  with  a  stronger  shock. 
The  vorticity  effects  are  greater,  and  as  before,  the 
potential  shock  is  downstream  of  the  zonal  solution, 
as  shown  in  Fig.  11.  The  Macb  contours  shown  in 
Figs.  12  and  13  correspond  to  the  potential  and  zonal 
solutions  respectively.  The  potential  shock  at  the  air¬ 
foil  trailing  edge  is  more  oblique  than  the  zonal  shock. 
The  velocity  defect  at  the  downstream  boundary  (x  ■= 
4)  is  given  in  Fig.  14  and  demonstrates  the  ability  of 
the  equations  to  convect  vorticity. 

As  a  final  example,  the  dual  potential  equations 
have  been  solved  in  general  curvilinear  coordinates 
{  (  =»  ((?■!/)■’!  ”  o(x,  y) )  for  a  subcritical  flow  over 
a  NACA  0012  airfoil  at  two  degrees  angle  of  attack. 
A  76  x  34  O-grid  was  used,  a  portion  of  which  is 
shown  in  the  vicinity  of  the  airfoil  (Fig.  15).  The  far 
field  is  rectangular  in  shape  with  rouoded  corners.  The 
outer  boundary  is  6  chord  lengths  from  the  airfoil  in 
the  x-direction  and  6  chord  lengths  in  the  y-direction; 
nevertheless  the  uniform  flow  boundary  condition 


i>  «  0  is  imposed  in  the  far  field.  The  value  of  ) 
on  the  body  is  obtained  from  the  Kutta  condition, 
Eq.  (32),  using  the  second-order  accurate  option.  A 
Cp  comparison  with  the  potential  solution  TAER7,t  is 
given  in  Fig.  16.  The  solutions  compare  very  well  even 
though  the  dual  potential  solution  is  computed  on  a 
coarser  grid  (50%  as  many  points).  The  convergence 
histories  of  the  ^  and  ^-residuals  are  shown  in  Fig. 
17. 

VL  CONCLUSIONS 

A  zonal  algorithm  which  utilizes  the  conservation- 
law-form  equations  together  with  the  dual  potential 
equations  has  been  developed  to  numerically  simulate 
steady  transonic  flow  of  an  inviscid  fluid  with  vor* 
ticity.  The  dual  potential  equations  are  able  to  con- 
vect  vorticity,  and  a  consistant  set  of  boundary  con¬ 
ditions  have  been  developed  which  permit  vorticity  at 
the  outflow  boundary.  The  zonal  algorithm  is  able  to 
capture  the  correct  shock  position.  The  dual  potential 
equations  require  significantly  less  CPU  time  than  the 
Euler  equations  since  they  are  solved  as  an  iteratively 
decoupled  set  of  scalar  equations.  The  test  problems 
we  have  investigated  showed  up  to  60%  savings  in  total 
CPU  time  for  a  zonal  solution  over  a  semi-flux  split ' 
solution  of  the  Euler  equations.  Code  optimization 
should  provide  additional  savings. 

The  dual  potential  equations  provide  an  alter¬ 
nate  method  of  computing  lifting  airfoil  flows.  Specifying 
^  on  the  body  in  order  to  satisfy  the  Kutta  condition 
eliminates  the  need  for  circulation  cuts;  however,  a 
Poisson  equation  must  be  solved.  The  increase  in  com¬ 
putational  work  is  perhaps  offset  by  coding  simplicity. 
This  will  be  especially  true  in  a  multi-element  prob¬ 
lem. 
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a)  Small  interactive  zone  about  a  shock 
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b)  Large  stationary  zone  about  a  shock 
Fig.  1  Zonal  partitioning. 
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Fig.  2  Dual  potential  boundary  conditions  for  asym 
metric  airfoil. 


Fig.  3  Zonal  boundary  condition  interfacing. 
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Fig.  6  Cp  comparison  between  semi-  flux  splitting  and 
centrally  differenced  Lulrr  scheme. 


Fig.  4  Kutta  condition  with  an  O-grid  topology. 
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Fig.  7  Cp  comparison  between  zonal  and  centrally 
differenced  Euler  scheme. 


^ig.  8  Cp  comparison  between  zonal  and  scalar  poten¬ 
tial  scheme. 
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Fig.  9  Convergence  history  of  scalar  potential  solu¬ 
tion. 
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Fig.  10  Convergence  histories  of  zonal  solution. 
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Fig.  15  O-grid  topology  near  NACA  0U12  airfoil. 


Fig.  16  Lifting  subcritical  Cp  comparison  between 
dual  potential  and  scalar  potential  scheme. 


Fig.  17  Convergence  histories  of  dual  potential  sola* 
tion. 
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